Ordinal compositional data and time series

3.1 Properties of ordinal distributions

Let Q be an ordinal RV with range S and PMF p ∈ S . To account for the natural order among the categories, one commonly does not focus on Q’s PMF but on its CDF, which can be summarized in the d -dimensional vector f = f 0 , … , f d - 1 ⊤ ∈ [ 0 ; 1 ] d . Here, we define f j = P Q ≤ s j for = 0 , … , d , where f d is omitted in f as it necessarily equals one. For this reason, we also used the largest category s d for ratioing in (2.1) when defining the alr transform.

Remark 2. A popular approach for a parametric modelling of ordinal data is the latent-variable approach (Agresti, 2010, p. 11). If f is a d -dimensional CDF vector, and if L is a (latent) realvalued R V with specified C D F F L ( x ) , then one represents f by selecting the threshold parameters - ∞ < η 0 < … < η d - 1 < + ∞ such that

f j = F L η j for j = 0 , … , d − 1

(3.1)

holds. In a sense, f is reparametrized in terms of η = η 0 , … , η d - 1 ⊤ . The most popular choice for F L is the (standard) logistic distribution, L ∼ L ( 0,1 ) , leading to the cumulative logit model. Here, L ( α , β ) denotes the logistic distribution with location parameter α and scale parameter β (see Johnson et al., 1995, Chapter 23), which is defined by the CDF

F α , β x = 1 + exp − x − α β − 1 for x ∈ ℝ .

(3.2)

The CDF f can be used to compute important stochastic properties of Q . The location of Q is commonly expressed in terms of the median (as this uses the order information), but also the mode (which is already well-defined for a nominal scale) might be used. Both median and mode take values in S , where the mode does not need to be unique. As S is typically of small size, the location measures provide only roughly graduated information about Q . Much more measures are available to quantify the dispersion of Q , see Blair and Lacy (1996), Kiesl (2003), Kvålseth (2011), Klein et al. (2016), Weiß 2019, 2022) for approaches and references. These measures are defined with respect to two extreme cases. Minimal dispersion happens if we are concerned with a one-point distribution on S , that is, if all probability mass concentrates in one category (maximal consensus):

f o n e ∈ t 0 , … , t d : = 1 1 ⋮ 1 , 0 1 ⋮ 1 , … , 0 ⋮ 0 1 , 0 ⋮ 0 0 ⊂ [ 0 ; 1 ] d.

By contrast, Q exhibits maximal dispersion if it follows the extreme two-point distribution (polarized distribution), f t w o = 1 2 , … , 1 2 ⊤ , that is, if we have 50 % probability mass in each of the outermost categories (maximal dissent). Note that the latter extreme differs from maximal nominal dispersion, which manifests itself in terms of a uniform distribution (maximal uncertainty). Ordinal dispersion measures are typically normalized to the range [ 0 ; 1 ] , and they take the value 0 (1) iff Q exhibits minimal (maximal) dispersion. A few examples (see the above references for many further options) are the index of ordinal variation,

I O V ( f ) = 4 d ∑ i = 0 d - 1 f i 1 - f i,

(3.3)

the coefficient of ordinal variation,

cov f = 1 - 1 d ∑ i = 0 d - 1 2 f i - 1 2 1 2 ;

(3.4)

and the cumulative paired entropy (using the convention 0 ⋅ l n 0 = 0 ),

C P E ( f ) = - 1 d l n 2 ∑ i = 0 d - 1 f i l n f i + 1 - f i l n 1 - f i.

(3.5)

Finally, the (a)symmetry of Q ’s distribution can be evaluated by measures of ordinal skewness, see Klein and Doll (2021), Weiß (2020) for details. We speak of symmetry if p i = p d - i for all i = 0 , … , d , or equivalently, if f i = 1 - f d - 1 - i for i = 0 , … , d - 1 . By contrast, maximal positive (negative) skewness is present if p 0 = 1 p d = 1 , that is, the extreme left (right) one-point distribution. This corresponds to the CDFs f = 1 = t 0 or f = 0 = t d , respectively, where 1 ( 0 ) denotes the vector of ones (zeros). These definitions of extreme skewness go along with the well-known rule of thumb that a positively skewed (= right-skewed) distribution appears to be left-leaning, and vice versa. A simple example of a normalized ordinal skewness measure is

skew f = 2 d ∑ i = 0 d − 1 f i − 1 with range − 1 ; 1 ,

(3.6)

OPEN IN VIEWER Figure 1

Age proportions in 195 countries from Example 1: plot of skew against IOV in (a), ternary diagrams in (b) and (c).

where skew = + 1 ( - 1 ) for maximal positive (negative) skewness, and where skew = 0 for a symmetric distribution.

To apply the ordinal measures of Section 3.1 to ordinal CoDa x ∈ S , we have to accumulate these vectors first, that is, to compute c : = T x with the matrix

T : = t 0 , … , t d = 1 0 ⋯ 0 ⋮ ⋱ ⋱ ⋮ 1 ⋯ 1 0 .

(3.7)

Then, we can evaluate the ordinal dispersion and skewness of x by applying formulae (3.3)-(3.6) to c instead of f . Let us illustrate this approach and its benefits by two real-world data examples.

Example 1. The dataset ageCat World from the R-package robCompositions contains the proportions x i = x i , 0 , x i , 1 , x i , 2 ⊤ , i = 1 , … , n , of people with age < 15 x i , 0 , 15 - 60 x i , 1 , and > 60 x i , 2 in n = 195 countries. These 3-part compositions refer to ordered categories, recall Section 1, thus constituting an example of ordinal CoDa. Since the dimension of the x i is three, a visual representation of the data by a ternary diagram is possible, see Figure 1. This shall help us to get an intuitive understanding of the ordinal measure (3.3)-(3.6) if being applied to the vectors c i : = T x i .

As indicated in Section 3.1, the location measures do not provide much discriminative power: in 191 out of 195 cases, the median category is s 1 (‘age 15-60’), and s 1 is also the mode category for 163 of the x i . Next, let us turn to ordinal dispersion. Here, we solely focus on the I O V c i values, as the IOV is most popular among practitioners, and as C O V c i and C P E c i provide essentially equivalent insights in the present example, in the sense of strong positive correlation to I O V c i , being > 0.99 in both cases (also see the discussion in Supplement S.1). As can be seen from Figure 1 (a), the values of I O V c i vary between 0.4 and 0.8, that is, we are concerned with medium to strong dispersion. In fact, only few countries lead to an IOV less than 0.6. These information are related to a ternary diagram of the data in Figure 1 (a). The few countries with IOV less than 0.6 (shown as empty circles) are those with a high proportion of people with age 15-60, that is, these compositions are quite close to a onepoint distribution (the three possible one-point distributions correspond to the corners of the triangle). Countries having an IOV ≥ 0.7 , in turn, are plotted (by grey dots) relatively close to the axis between the lowest and largest age category, that is, the outermost categories are more dominant.

Next, we consider the s k e w c i values. Figure 1 (a) shows that we have at most moderate skewness values, with maximal absolute extent around 0.4, where positive skewness values are slightly more frequent than negative ones. Note that countries with negative skewness tend to show larger dispersion than in the case of positive skewness. The ternary diagram in Figure 1 (c) shows that positive skewness values refer to the lower-left half of the triangle, which is easy to interpret: right-skewed compositions are left leaning, so the lower categories dominate the upper ones.

Example 1 showed that for ordinal 3-part compositions, the IOV-skew diagram provides similar insights into the data as the ternary diagram. However, an IOV-skew diagram can also be drawn for higher-dimensional compositions, whereas a ternary diagram is not immediately applicable (one could combine individual categories into larger ones before plotting). By contrast, the IOV-skew diagram allows to directly visualize the high-dimensional ordinal compositions, being easy to interpret by a practitioner, see the data example in Supplement S.1. To sum up, Section 3 demonstrated that ordinal measures such as IOV and skew, if applied to ordinal CoDa, provide new insights into the data (which may also be presented visually) and, thus, constitute a valuable complement to existing methods for descriptive CoDa analysis.

4.1 Asymptotic distributions and simulation results

Let X 1 , … , X n be an independent and identically distributed (i.i. d.) sample of ordinal compositional RVs with mean p , let X ‾ denote their sample mean. Then, the central limit theorem (CLT) implies that n ( X - - p ) is asymptotically normally distributed (see Serfling, 1980), namely

n X ‾ − p ∼ N 0 , Σ with Σ = σ i j i , j = 0 , … , d and σ i j = Cov X i , X j .

(4.1)

While (4.1) holds for nominal CoDa as well, we now turn to the accumulated versions of X and p , that is, to C k = T X k for k = 1 , … , n with mean E C k = f = T p , where T is the transformation matrix from (3.7). These are meaningful only for ordinal CoDa, and they are required for the computation of the ordinal statistics (3.3)-(3.6). The asymptotics of the sample mean C - follow immediately from (4.1) and the linearity properties of the multivariate normal distribution, namely

n ( C - - f ) ∼ N 0 , Σ ' w i t h Σ ' = T Σ T ⊤ .

(4.2)

The entries of Σ ' are computed as

σ i j ′ = T i • Σ T j • ⊤ = ∑ r = 0 i ∑ s = 0 j σ r s for i , j = 0 , … , d − 1 ,

where T k ⦁ abbreviates the k th row vector of T . Finally, the asymptotics of the ordinal statistics (3.3)(3.6), if applied to the estimator C ¯ , are derived based on Taylor expansions (‘Delta method’), see Supplement S.2 for details. The obtained results are summarized in the following theorem.

Theorem 1. Let C ¯ be computed from an i.i.d. ordinal CoDa sample, let the σ i j ′ be given as in (4.2). Then, the ordinal measures (3.3)-(3.6) are asymptotically normally distributed with asymptotic variances

V [ IOV ( C ¯ ) ] ≈ 1 n 16 d 2 ∑ i , j = 0 d − 1 1 − 2 f i 1 − 2 f j σ i j ′

V [ COV ( C ¯ ) ] ≈ 1 n 4 d 2 ∑ i , j = 0 d − 1 1 − 2 f i 1 − 2 f j ( 1 − COV ) 2 σ i j ′ ;

V [ CPE ( C ¯ ) ] ≈ 1 n 1 d 2 ( ln 2 ) 2 ∑ i , j = 0 d − 1 ln 1 − f i f i ln 1 − f j f j σ i j ′ ;

V [ skew ( C ¯ ) ] ≈ 1 n 4 d 2 ∑ i , j = 0 d − 1 σ i j ′ .

While E [ s k e w ( C - ) ] = s k e w ( f ) holds exactly, the means of the remaining measures are approximately equal to

E [ IOV ( C ¯ ) ] ≈ IOV ( f ) − 1 n 4 d ∑ i = 0 d − 1 σ i i ′ ;

E [ COV ( C ¯ ) ] ≈ COV ( f ) + 1 n 2 d ∑ i , j = 0 d − 1 1 − 2 f i 1 − 2 f j d ( 1 − COV ) 3 − δ i j 1 − C O V σ i j ′ ;

E [ CPE ( C ¯ ) ] ≈ CPE ( f ) − 1 n 1 2 d ln 2 ∑ i = 0 d − 1 1 f i 1 − f i σ i i ′ .

OPEN IN VIEWER Figure 2

Proportions of education levels in 31 European countries from Section 4.2: ternary diagram in (a), plot of skew against IOV in (b), estimates in (c).

For the practical implementation of the asymptotics, it is suggested to compute the σ i j ' separately, and to plug-in the results into the formulae of Theorem 1. In the special case of a Dirichlet distribution, see Example 1 in Supplement S.2, closed-form formulae for the σ i j ' exist such that further simplifications are possible.

In practice, the asymptotic distributions of Theorem 1 are used to approximate the true distributions of the ordinal statistics IOV C ¯ , … , skew C ¯ according to (3.3)-(3.6). Therefore, it is important to check the approximation quality of these distributions for finite sample sizes n . For this purpose, the simulation study summarized in Supplement S.3 was done, where i. i. d. ordinal CoDa samples according to a Dir p , v -distribution and a LN μ , Σ -distribution were simulated (with 10 6 replications per scenario). The actual parametrization of these distributions is inspired by the subsequent data example in Section 4.2. Details on the simulation procedures as well as on required numerical approaches for the LN-case are discussed in Supplement S.3. It turns out that we generally have an excellent agreement between simulated and asymptotic properties. Thus, the bias and SD being observed in finite samples are well described by the asymptotics of Theorem 1 .

We analyse the dataset educFM in the R-package robCompositions (also see Filzmoser et al., 2018, Section 10.7.3), which provides the proportions of low, medium, and high education levels of father and mother in n = 31 European countries. So for each country i ∈ 1 , … , n , two 3-part compositions are available, say x f , i and x m ; i for father and mother, respectively. The category s 0 (‘low’) is the median (mode) of x f ; i in 19 (22) cases, and the median (mode) of x m ; i in 22 (25) cases; otherwise, it is the category s 1 (‘medium’). The ternary plot in Figure 2 (a) shows that the CoDa concentrate in the low-to-medium corner, which goes along with the positive skew c * ; i values in part (b), implying left-leaning distributions. The IOV c * ; i values cover a wide range from low to high dispersion. It also interesting to note the pronounced structure in the IOV-skew diagram, where low dispersion is connected to strong positive skewness, that is, the corresponding CoDa vectors are close to a one-point distribution in the lowest category.

Next, we analyse if there are significant differences (on a 5%-level) between the education proportions of father and mother; from the plots in Figure 2 (black vs. grey dots), no such difference can be recognized. For this purpose, models for both datasets are required. A natural first idea is to try a Dirichlet distribution, but a problem becomes clear: the sample covariance matrices for the fathers’ and mothers’ compositions are 0.036 − 0.027 − 0.009 − 0.027 0.022 0.005 − 0.009 0.005 0.005 and 0.029 − 0.022 − 0.007 − 0.022 0.019 0.003 − 0.007 0.003 0.004 , respectively, that is, some cross-covariances are positive, which is not consistent with a Dirichlet distribution, recall (2.4). Thus, logistic normal distributions seem more appropriate, as also confirmed by the distributional analyses presented in Supplement S.4. Thus, the following discussion focusses on the LN model (fitted by mvnorm.mle in the R-package Rfast), whereas results regarding the Dirichlet model are presented in Supplement S.4.

The fitted LN μ , Σ -models are μ f = ( 1.451 , 0.797 ) ⊤ , Σ f = 0.721 0.028 0.028 0.278 and μ m = ( 1.969 , 1.084 ) ⊤ , Σ m = 0.728 0.171 0.171 0.418 , respectively, which both lead to a rather good fit to the CoDa's means and cross-covariances. Hence, the LN-models are used as the base for asymptotic computations. The point estimates from IOV C ¯ and skew C ¯ , as well as the corresponding asymptotic biases and standard errors (SEs) according to Theorem 1, are shown in Figure 2 (c). The biases of IOV C ¯ are negligible relative to the actual estimates, and skew is unbiased anyway. Comparing the IOV and skew estimates between father and mother, we recognize a notable difference, which is larger than twice a standard error. So in view of the estimators’ asymptotic normality, we conclude that the IOV and skew estimates are significantly different between father and mother. In fact, we have less dispersion and stronger positive skewness for mothers, that is, the mothers’ education proportions are (in the mean) closer to a one-point distribution in the category ‘low’. This might be interpreted as an indication of unequal opportunities for education among males and females in the past.

5.1 Control charts for ordinal CoDa

In what follows, we are concerned with the monitoring of an ordinal CoDa process X t t ∈ ℕ , which is assumed to be i. i. d. in time t under in-control conditions. There are several proposals for CoDa control charts in general (i. e., nominal CoDa), see Holmes and Mergen (1993), Vives-Mestres et al. 2013, 2014 and references therein, but none for the particular case of ordinal CoDa. In fact, the data application (manufacturing of grit) considered in the aforementioned articles exhibits ordered categories, because it refers to the particle size proportions with categories ‘large’, ‘medium’, and ‘small’; see Table 2 in Holmes and Mergen (1993) for a sequence of 56 CoDa vectors. This ordinal information, however, was not used for defining the control charts. Instead, Holmes and Mergen (1993), Vives-Mestres et al. 2013, 2014 used types of Hotelling- T 2 chart (see Montgomery, 2009, Section 11.3), which apply to unordered categories as well: while Holmes and Mergen (1993) computed the T 2 statistic based on two of the three components of X t , Vives-Mestres et al. (2013, 2014) first mapped the 3-part simplex S onto ℝ 2 and then computed the T 2 -statistic. Generally, see Section 11.3 in Montgomery (2009), if Y t is a multivariate vector-valued process with in-control mean μ 0 and covariance matrix Σ 0 , the χ 2 -chart is defined by plotting the statistics

X t 2 = Y t − μ 0 ⊤ Σ 0 − 1 Y t − μ 0 for t = 1 , 2 , …

(5.1)

against an upper CL (UCL) chosen such that the ARL 0 -target is met, that is, an alarm is triggered if X t 2 > UCL . The Hotelling- T 2 chart differs from (5.1) by using estimates for μ 0 and Σ 0 instead (and by adapting the UCL appropriately). Besides not using the ordinal information, statistic (5.1) has the drawback that an alarm is difficult to interpret (see Montgomery, 2009, p. 505), because the value of X t 2 does indicate the actual type of out-of-control situation; see Vives-Mestres et al. (2013) for a discussion. Thus, if monitoring ordinal CoDa, it is suggested to complement the CoDa χ 2 - or T 2 -chart by further control charts based on the ordinal statistics in Section 3.1, such that a targeted diagnosis of the actual out-of-control state is possible.

To monitor the process X t or C t , respectively, two types of IOV-chart are developed by adapting the corresponding proposals in Ottenstreuer et al. (2023). Let p 0 and f 0 denote the in-control means of X t and C t , respectively. We shall consider two-sided chart designs, that is, with both a lower and upper CL (LCL and UCL, respectively). Then, the Shewhart IOV-chart is defined by plotting

IOV C t for t = 1 , 2 , …

(5.2)

against LCL and UCL, which are chosen such that the in-control ARL is close to A R L 0 . An alarm is triggered if either IOV C t < LCL or IOV C t > UCL . As this decision relies solely on the t th observation C t , Shewhart-type charts typically react quickly to very large process shifts, but rather slow to small process changes.

To equip the IOV-chart with an additional memory and to, thus, improve the sensitivity towards small process changes, we also consider an exponentially weighted moving-average (EWMA) approach in analogy to Ottenstreuer et al. (2023). For the smoothing parameter λ ∈ 0 ; 1 , we recursively compute C t , λ = λ C t + 1 − λ C t − 1 , λ with C 0 , λ = f 0 . Then, the EWMA IOV-chart is defined by plotting

IOV C t , λ for t = 1 , 2 , …

(5.3)

against appropriately chosen LCL, UCL. Note that λ → 1 leads to the Shewhart chart (5.2), whereas decreasing values of λ increase the inherent memory of (5.3).

Both IOV-charts (5.2) and (5.3) are defined to detect changes in the ordinal dispersion. A violation of LCL (UCL) indicates a decrease (increase) of dispersion and, thus, allows for a targeted diagnosis. Instead of the IOV (3.3), one could also use the COV (3.4) or CPE (3.5), but we focus on the IOV for convenience. However, it seems reasonable to also consider the Shewhart and EWMA version of a skew chart based on (3.6), as such a chart could react to changes in ordinal skewness, see Ottenstreuer et al. (2023). Hence, we also plot

skew C t , λ for t = 1 , 2 , …

(5.4)

against appropriate LCL , UCL , where λ = 1 corresponds to the Shewhart skew-chart, and λ ∈ 0 ; 1 to the EWMA skew-chart. To sum up, at least theoretically, the IOV- and skew-charts offer the potential for a targeted diagnosis of the actual out-of-control scenario. In the subsequent Section 5.2, we provide empirical evidence that such targeted conclusions are indeed possible.

As mentioned in Section 5.1, Vives-Mestres et al. 2013, 2014 applied a T 2 -type control chart to the CoDa process X t , namely after first mapping the 3-part simplex S onto ℝ 2 . A possible solution for this mapping could have been the alr transform 2.1 , but Vives-Mestres et al. 2013, 2014 used types of ilr transforms instead. To make the subsequent results comparable to those in Vives-Mestres et al. (2014), in this section, we shall consider their specific ilr transform defined by

ilr x = y ∈ ℝ 2 with y 1 = 1 2 ln x 1 x 2 , y 2 = 2 3 ln x 1 x 2 x 0 .

(5.5)

It should be noted, however, that because of the relation ln x 1 x 2 x 0 = 1 2 ln x 1 x 2 − ln x 0 x 2 , it holds that ilr x = 0 1 / 2 − 2 / 3 1 / 6 alr x , that is, ilr and alr only differ by a linear transformation, also recall Remark 1. We use the CoDa-normal distribution implied by requiring that ilr X ∼ N μ , Σ with μ ∈ ℝ 2 and positive definite Σ ∈ ℝ 2 × 2 , also see Mateu-Figueras et al. (2013), which is a reparametrization of the LN-distribution discussed in Section 2. This distribution is relevant for chart design, as we have to specify the in-control model for the process X t . Like in Vives-Mestres et al. (2014), we define the in-control model as ilr X ∼ N μ 0 , Σ 0 with μ 0 ≈ ( 1.959 , 1.158 ) ⊤ and Σ 0 ≈ 0.134 − 0.032 − 0.032 0.097 , which are the sample values from the ilr-transformed particles dataset x 1 , … , x 56 in Table 2 of Holmes and Mergen (1993). For comparability, we also use the same target ARL , namely ARL 0 = 333. 3 ¯ . Note that the above μ 0 corresponds to the component-wise means p 0 ≈ ( 0.057 , 0.881 , 0.062 ) ⊤ (computed by numerical integration) and component-wise medians ilr^-1 μ 0 ≈ ( 0.054 , 0.890 , 0.056 ) ⊤ for X , that is, we have a nearly symmetric mean/median particle size proportion under in-control conditions, which also exhibits only low dispersion. More precisely, IOV p 0 ≈ 0.223 and skew p 0 ≈ − 0.005 . Note that a monitoring of ordinal location would not be promising: because of the large mean/median proportion of category s 1 (‘medium’), we can hardly ever expect to observe any other median or mode category than s 1 . In fact, for the data in Table 2 of Holmes and Mergen (1993), s 1 is the median and mode category without exception.

The aforementioned in-control model is now used for chart design. For given values of the design parameters, 10 5 in-control processes are simulated and the run lengths (i. e., the numbers of statistics until the first signal) are determined; the mean across these run lengths approximates the true ARL value. For chart design, the design parameters are varied (and ARL simulation is repeated for each revised design) until the ARL 0 -target is met in close approximation. In this way, we determine the chart design of the Shewhart and EWMA versions of IOV- and skew-chart (5.2)-(5.4), where the EWMA charts use λ = 0.1 . As a competitor, we consider the Shewhart χ 2 -chart (5.1) in analogy to Vives-Mestres et al. (2014). Furthermore, for the sake of completeness, we also consider the multivariate EWMA version of this chart defined like in Section 11.4 of Montgomery (2009):

X t , λ 2 = 2 − λ λ 1 − ( 1 − λ ) 2 t − 1 Z t , λ − μ 0 ⊤ Σ 0 − 1 Z t , λ − μ 0 for t = 1 , 2 , … ,

(5.6)

where Z t , λ = λ Y t + 1 − λ Z t − 1 , λ with Z 0 , λ = μ 0 . Although both EWMA charts (5.3) and (5.6) are not directly comparable because of a different implementation of the EWMA approach, we use again λ = 0.1 for illustration. The resulting designs and in-control ARLs are summarized in the first block (italic font) of Table S.6.1 in Supplement S.6.

Next, these chart designs are used for a simulation-based analysis of the out-of-control performance. Based on an ordinal logit approach (recall Remark 2), we define several out-of-control scenarios for ilr^-1 (μ) (medians of CoDa distribution), see Supplement S.5 for details. The logit approach accounts for the order among the categories, and it allows to distinguish between dispersion and skewness shifts for ilr − 1 μ compared to ilr − 1 μ 0 . The simulation results discussed in Supplement S. 6 show that the IOV-, skew-, and χ 2 -charts have a quite different out-of-control performance. While the χ 2 -chart is universally applicable but does not allow to conclude from the signal on the type of out-of-control scenario, the IOV-and skew-charts are sensitive only to specific out-of-control scenarios, namely such affecting dispersion and skewness, respectively. Thus, running an IOV- and skew-chart in parallel to a χ 2 -chart (preferably always the EWMA version), we can expect that a dispersion shift leads to an early alarm of the IOV-chart, and a skewness shift to one of the skew-chart. In this way, the types of alarm allow us to conclude on the type of out-of-control situation.

This type of diagnostic feature is illustrated by Figure 3, where the EWMA IOV-and skew-chart are applied to the particle size proportions of Table 2 in Holmes and Mergen (1993). The skewchart in (b) does not trigger any alarm during the monitoring period, so we have no indication of a systematic skewness change. The IOV-chart in (a), by contrast, triggers alarms at times t = 16 , 17 as well as t = 23 , 24 . In both cases, the LCL is violated, indicating that there was an exceptionally low dispersion for t ≤ 16 and t ≤ 23 . Further analysing the raw data, we indeed note that x 15 = 0.026 , 0.945 , 0.029 and x 16 = 0.027 , 0.945 , 0.028 are considerably less dispersed than ilr − 1 μ 0 ≈ ( 0.054 , 0.890 , 0.056 ) ⊤ (but still nearly symmetric). Note that x 15 , x 16 are close to the second scenario in the fifth block of Table S.6.1. Analogous results (though less pronounced) hold for observations x 21 , … , x 24 . Finally, although it did not suffice for an alarm yet, we recognize increasing dispersion and skewness towards the end of the monitoring period in Figure 3.

OPEN IN VIEWER Figure 3

Particle size proportions from Section 5.2: EWMA IOV-chart in (a) and EWMA skew-chart in (b), where λ = 0.1 .

6.1 Analysis and modelling of ordinal CoTS

Several models for CoDa processes X t and the resulting CoTS X 1 , … , X n have been proposed during the last years. The most common approach for modelling CoTS is to first map the d + 1 part compositions X t onto d -dimensional real vectors Y t (e. g., by taking the alr transform (2.1) or an ilr transform), and to then apply a vector autoregressive moving-average (VARMA) model to the transformed time series Y t , see, for example, Mills (2010), Dawson et al. (2014). Because of the large number of VARMA model parameters, for short time series, also simple regression approaches have been proposed, that is, the component series Y t , i of the transformed CoTS are considered as linear or quadratic regression models in time t , see, for example, Wang et al. (2007), Mills (2010). A parsimonious ARMA-type model for CoTS relying on the use of a variation operator was proposed by Möller and Weiß (2020). Also Zheng and Chen (2017) proposed a class of ARMA-type models for CoTS, being based on a conditional regression approach. For references on further model types, see Pawlowsky-Glahn and Buccianti (2011, Chapter 7) and Zheng and Chen (2017).

All the aforementioned time series models are developed to deal with general (nominal) CoDa, that is, a possible order among the categories is not considered. However, several of the examples presented in the literature are ordinal CoTS. For example, Wang et al. (2007) discuss a time series of yearly proportions of employed persons in China, where the proportions refer to the three production sectors ‘primary industry’, ‘secondary industry’, and ‘tertiary industry’; also see Supplement S.7. Mills (2010) modelled obesity trends in England, expressed by yearly proportions of persons falling into ordered weight categories, ranging from ‘not overweight’ to ‘obese’; also see Section 6.2. As a last example, Dawson et al. (2014) analysed a time series expressing the attitude of UK people towards the Olympic Games, with ordered categories from ‘strongly against’ to ‘strongly supportive’. In any of these examples, it could be beneficial to explicitly account for the ordered categories while analysing the data, which is illustrated by the data example discussed in Supplement S.7.

Not only for their analysis, but also for the modelling of ordinal CoTS, it could be relevant to account for the ordered categories. As a possible solution, we combine the conditional regression model of Zheng and Chen (2017) with the cumulative logit approach of Remark 2, in analogy to the ordinal time series models of Pruscha (1993); Fokianos and Kedem (2003). Zheng and Chen (2017) propose a conditional modelling of X t . If F s denotes the σ -field generated by all information up to time s , then X t ∣ F t − 1 is assumed to follow a specified CoDa distribution with parameter vector θ t , say X t ∣ F t − 1 ∼ CoDa θ t , where θ t is computed from the information up to time t − 1 . For example, we would define θ t = p t , 0 , … , p t , d − 1 , v ⊤ if X t is conditionally Dirichlet-distributed with time-invariant scale parameter v , that is, X t ∣ F t − 1 ∼ Dir p t , v according to (2.3), where p t , d = 1 − p t , 0 − … − p t , d − 1 . In this case, the conditional mean of X t equals E X t ∣ F t − 1 = p t , and the conditional (co)variances are adjusted by the factor 1 / 1 + v , recall (2.4). As suggested by Zheng and Chen (2017), p t could be modelled by an ARMA-type recursion of the form p t = ϕ X t − 1 , … , X t − p , p t − 1 , … , p t − q , where p ∈ ℕ expresses the autoregressive order and q ∈ ℕ 0 = 0 , 1 , … the ‘moving-average’ order (it is more appropriate to interpret p t − 1 , … , p t − q as forming a feedback term). Specific examples for choosing ϕ in the case of nominal CoDa are discussed by Zheng and Chen (2017).

In what follows, to account for the natural order of ordinal CoDa and to reduce the number of model parameters at the same time, we present possible specifications for ϕ being derived from the cumulative logit approach, in analogy to Pruscha (1993); Fokianos and Kedem (2003). Define X ˜ t = X t , 0 , … , X t , d − 1 ⊤ and p ˜ t = p t , 0 , … , p t , d − 1 ⊤ , and let f t denote the corresponding d -dimensional conditional CDF vector. Recall that F 0 , 1 denotes the CDF of the standard logistic distribution (inverse logit function). An ARMA-type model for ordinal CoTS is defined by the recursion

f t , i = F 0 , 1 η i + ∑ k = 1 p α k ⊤ X ˜ t − k + ∑ l = 1 q β l ⊤ p ˜ t − l for i = 0 , … , d − 1 ,

(6.1)

which has the d p + q + 1 parameters involved in η , α 1 , … , α p , β 1 , … , β q . The parameters α 1 , … , β q control the serial dependence structure, which is unique across the categories. The individual ordinal categories are accounted for by the parameter η = η 0 , … , η d − 1 ⊤ with − ∞ < η 0 < … < η d − 1 < + ∞ . A more parsimonious model specification, also see Pruscha (1993), is

f t , i = F 0 , 1 η i + ∑ k = 1 p α k 1 ⊤ X ˜ t − k + ∑ l = 1 q β l 1 ⊤ p ˜ t − l for i = 0 , … , d − 1 ,

(6.2)

which has only d + p + q parameters. If the data are then generated according to X t ∣ F t − 1 ∼ Dir p t , v , for example (where p t follows from f t by taking discrete differences), we also have the additional scale parameter v from the Dirichlet distribution. Models (6.1) and (6.2) are easily adapted to account for possible covariate information z t by adding the summand ‘ + γ ⊤ z t ’ within the parentheses.

OPEN IN VIEWER Figure 4

Proportions of weight categories from Section 6.2: (a) plot of proportions over time, and (b) table with model fits and corresponding CSS values.

Parameter estimation can be done by a conditional least-squares (CLS) approach, that is, by numerically minimizing the conditional sum of squares (CSS) defined by

CSS ( θ ) : = ∑ t ∥ X t − p t ∥ 2 → min ,

(6.3)

where θ comprises all parameters involved in (6.1) or (6.2), respectively, and where ∥ ⋅ ∥ denotes the Euclidean norm. The CLS approach is semi-parametric in the sense that it does not require to specify the full conditional distribution of X t ∣ F t − 1 but only the equation for the conditional mean, which, in turn, is useful for predicting future CoDa values. It is also easy to formulate a conditional log-likelihood approach for models (6.1) and (6.2), but this requires to fully specify the distribution of X t ∣ F t − 1 and is computationally more demanding (due to the complex structure of typical CoDa distributions).

The following data example is inspired by the one in Section 4.2 of Mills (2010), who modelled obesity trends in England for the years in 1993-2005. We consider, however, a much longer period, namely yearly proportions for the period 1975-2016. More precisely, based on the (agestandardized) percentages for different BMI classes (body mass index) offered by the World Health Organization at https://apps.who.int/gho/data/node.main.BMIANTHROPOMETRY?lang=en, the CoTS x 1 , … , x 42 for Germany is derived, with the three categories (hence d = 2 ) ‘not overweight’ (BMI < 25), ‘overweight’ (BMI in [25;30)), and ‘obese’ (BMI ≥ 30 ). The data are plotted in Figure 4 (a), where the proportions correspond to the heights of the grey areas (the curves between the grey areas refer to the cumulative proportions). The vertical dashed line between the years 2010 and 2011 expresses the following partition: the data x 1 , … , x 36 for 1975-2010 are used for model fitting, whereas the data x 37 , … , x 42 for 2011-2016 are left for out-of-sample forecasting.

Looking at Figure 4 (a), an obvious first candidate for modelling the CoTS is a simple (logistic)linear model, that is, (6.1) with

f t , i = F 0 , 1 η i + γ i t for i = 0 , … , d − 1 .

(6.4)

This is analogous to the model used by Mills (2010) for their UK series, being linear in the components of the alr-transformed data:

E alr x t i = a i + b i t for i = 1 , … , d .

(6.5)

Model (6.5) is used as a competitor here. Note that (6.5) could be equivalently expressed in terms of the ilr transform (5.5), which would lead to a reparametrization of the model. Finally, we also consider a version of (6.1) having both a linear trend and a first-order AR-component, namely

f t , i = F 0 , 1 η i + γ i t + α 1 ⊤ X ˜ t − 1 for i = 0 , … , d − 1 .

(6.6)

The results of CLS estimation (based on x 1 , … , x 36 ) are summarized in Figure 4 (b). The linear coefficients for the fitted ordinal models (6.4) and (6.6) are negative, which is reasonable in view of the decreasing curves for c t , 0 = x t , 0 and c t , 1 = x t , 0 + x t , 1 in Figure 4a. The negative linear coefficients for model (6.5) refer to the alr transform (2.1) of the non-cumulative proportions, so they indicate decreasing quotients x t , 0 / x t , 2 , x t , 1 / x t , 2 (if we would use the ilr (5.5) instead of the alr, the parameter estimates of a i , b i would change to 0.869 , − 0.013 for i = 1 and − 1.096 , 0.019 for i = 2 . It is also interesting to note that the estimates for η , γ are quite similar between models (6.4) and (6.6), that is, one may interpret the estimate of the AR-coefficient α 1 in (6.6) as a refinement of the baseline linear model (6.4). To compare the performance of the fitted models, two types of CSS are computed: ‘ CSS in ’ equals 10 3 times the CSS (6.3) computed for the in-sample data t = 2 , … , 36 (the summand t = 1 is left out for all models for a fair comparison, because this summand is not defined for the AR model (6.6)), and ‘ CSS out ’ equals 10 3 times the CSS (6.3) computed for the out-of-sample data t = 37 , … , 42 . Here, the unique scaling factor 10 3 was used to improve the readability of the printed CSS values. It can be seen that the ordinal linear model (6.4) clearly outperforms the linear alrmodel (6.5), because both CSS values of (6.4) are equal to only about 60 % of the (6.5)’s CSS values (the CSS values of the ilr-version of model (6.5) are identical to those of the alr-version). However, a further improvement (especially in terms of the out-of-sample performance) is achieved by the additional AR-term in (6.6). This term, given by ‘ + 0.322 X t − 1 , 0 − 0.651 X t − 1 , 1 ’, shows that the past proportion for ‘not overweight’ has an increasing effect on f t , 0 , f t , 1 , whereas the one of ‘overweight’ has a decreasing effect. As x t , 0 strongly decreases for t = 1 , … , n while x t , 1 is rather stable, the AR-term becomes increasingly negative over time and, thus, intensifies the linearly decreasing trend in F 0 , 1 − 1 f t , i .