MUFF, STEFANIE, RIEBLER, ANDREA, HELD, LEONHARD, RUE, HAVARD ˚ , & SANER, PHILIPPE. 2015. Bayesian analysis of measurement error models using integrated nested Laplace approximations. *Journal of the*

SIMPSON, DANIEL, RUE, HAVARD ˚ , RIEBLER, ANDREA, MARTINS, THI-AGO G, & SØRBYE, SIGRUNN H. 2017. Penalising model component complexity: A principled, practical approach to constructing priors. *Sta-*

*Royal Statistical Society: Series C: Applied Statistics*, 231–252. RIEBLER, ANDREA, SØRBYE, SIGRUNN H, SIMPSON, DANIEL,&RUE, HAVARD ˚ . 2016. An intuitive Bayesian spatial model for disease mapping that accounts for scaling. *Statistical methods in medical research*, 25(4),

Figure 1. *Map displaying the posterior means of car crashes counts.*

β*x* 0.01 1.064 2.95 β0 -5.307 -9.90 -15.441 βprimary 0.61 0.56 0.40 βsecondary 0.57 0.68 1.05 DIC 33126 30466 Table 1. *Summary of DIC, posterior means of fixed effects, and error-prone covariate.*

statistical analysis.

References

1–101.

*arXiv:2011.12595*.

1145–1165.

*tistical science*, 1–28.

where β0 denotes the intercept, *zi* is an error free covariate representing the road-type of each segment, *xi* is an unobservable error prone covariate representing the traffic volumes, while β*x* and β*z* are the corresponding coefficients. Finally, θ*i* and φ*i* denote spatially structured and unstructured random effects that are modelled using a reparametrisation and a network re-adaptation of Besag-York-Mollie (BYM) prior (Riebler ´ *et al.*, 2016, Gilardi *et al.*, 2020).

The classical spatial ME model assumes that *xi* can be observed only via a proxy, say *wi*, such that

$$\boldsymbol{w}\_{i} = \boldsymbol{x}\_{i} + \boldsymbol{\mu}\_{i} + \boldsymbol{\Phi}\_{i}; \ i = 1, \ldots, n.$$

The terms *ui* and ϕ*i* represent the ME and denote, respectively, spatially structured and unstructured random effects that are also modelled using the BYM prior. In particular, parameter ϕ*i* adds a spatial smoothing effect to the unobserved covariate *xi*. In this paper, we assume that the edge betweenness centrality measure can approximate the unobservable traffic volumes.

At the second stage of the hierarchy, we specified an exposure model that relates *xi* with the error-free predictor:

$$\mathfrak{a}\_{i} = \mathfrak{a}\_{0} + \mathfrak{a}\_{\mathfrak{z}} \mathfrak{z}\_{i} + \mathfrak{e}\_{i};\ i = 1, \ldots, n. \tag{2}$$

The parameter α0 denotes the intercept, α*z* is the coefficient of the error-free covariate, and ε*i* is a normally distributed error component. Furthermore, we assigned independent *N*(0,103) priors to β0, β*z*, α0, and α*z*, i.e. the intercepts and the coefficients assigned to *zi* in equations (1) and (2).

The third level completes the specification of the hierarchical model eliciting a N(0,100) prior for β*x*, i.e. the coefficient of the error-prone covariate, a Gamma(1,5e-05) prior on the precision of *ui* and ϕ*i*, and Penalised Complexity priors for the parameters of BYM's re-adaptation (Simpson *et al.*, 2017).

#### 4 Results and conclusions

We estimated the statistical model described in Section 3 using INLA methodology and compared the results with two simpler models: the first one completely ignores ME, while the second one adopts a classical ME without spatial smoothing effects. We found that omitting ME greatly attenuates the importance of traffic volumes, and excluding the spatial smoothing terms worsens

Table 1. *Summary of DIC, posterior means of fixed effects, and error-prone covariate.*

Figure 1. *Map displaying the posterior means of car crashes counts.*

the fit of the model. Motorways were found less prone to car crashes than the other road types, while the posterior distributions of fixed effects and common hyperparameters were found stable among the three models. We report in Table 1 a short summary of fixed effects' posterior means, while Figure 1 displays the posterior means of predicted counts. We can notice that it highlights a few road segments close to the city centre that would require a more detailed statistical analysis.

#### References

In the first level of the hierarchy, we define a log-linear structure on λ*i*, i.e.

where β0 denotes the intercept, *zi* is an error free covariate representing the road-type of each segment, *xi* is an unobservable error prone covariate representing the traffic volumes, while β*x* and β*z* are the corresponding coefficients. Finally, θ*i* and φ*i* denote spatially structured and unstructured random effects that are modelled using a reparametrisation and a network re-adaptation of Besag-York-Mollie (BYM) prior (Riebler ´ *et al.*, 2016, Gilardi *et al.*, 2020). The classical spatial ME model assumes that *xi* can be observed only via a

*wi* = *xi* +*ui* +ϕ*i*; *i* = 1,...,*n*.

The terms *ui* and ϕ*i* represent the ME and denote, respectively, spatially structured and unstructured random effects that are also modelled using the BYM prior. In particular, parameter ϕ*i* adds a spatial smoothing effect to the unobserved covariate *xi*. In this paper, we assume that the edge betweenness

At the second stage of the hierarchy, we specified an exposure model that

The parameter α0 denotes the intercept, α*z* is the coefficient of the error-free covariate, and ε*i* is a normally distributed error component. Furthermore, we assigned independent *N*(0,103) priors to β0, β*z*, α0, and α*z*, i.e. the intercepts

The third level completes the specification of the hierarchical model eliciting a N(0,100) prior for β*x*, i.e. the coefficient of the error-prone covariate, a Gamma(1,5e-05) prior on the precision of *ui* and ϕ*i*, and Penalised Complexity priors for the parameters of BYM's re-adaptation (Simpson *et al.*, 2017).

We estimated the statistical model described in Section 3 using INLA methodology and compared the results with two simpler models: the first one completely ignores ME, while the second one adopts a classical ME without spatial smoothing effects. We found that omitting ME greatly attenuates the importance of traffic volumes, and excluding the spatial smoothing terms worsens

*xi* = α0 +α*zzi* +ε*i*; *i* = 1,...,*n*. (2)

centrality measure can approximate the unobservable traffic volumes.

and the coefficients assigned to *zi* in equations (1) and (2).

proxy, say *wi*, such that

relates *xi* with the error-free predictor:

4 Results and conclusions

log(λ*i*) = β0 +β*zzi* +β*xxi* +θ*i* +φ*i*; *i* = 1,...,*n*, (1)

## THE *Lp* DATA DEPTH AND ITS APPLICATION TO MULTIVARIATE PROCESS CONTROL CHARTS

Multivariate quality control studies was first conducted by Hotteling, 1947.

Woodall, 2000 distinguishes the techniques of the control chart processing in Phase I (also called retrospective or preliminary phase) and in Phase II. Phase I uses charts with the purpose of defining whether a process is statistically under control when the first group are processed. In Phase II, the charts are used to check if the process is in control when future subgroups were being processed. In this last phase, it is assumed that the distribution of the process is known and most of the classical applications require the hypothesis that the process under consideration follows a multivariate normal distribution. However, in most industrial applications the distribution of a parametric multivariate control chart is difficult to estimate for processes with multiple quality characteristics. As such, all observations are considered as *d*-dimensional vector and therefore used to detect possible shifts in the *d*-dimensional distributions of the quality process. A statistical process control (SPC) procedure set up in a multivariate framework is more effective than a joint monitoring system consisting of a se-

The most popular multivariate statistical process control charts are based on the Hotelling's *T*2 statistics, that are a multivariate extension of Shewart's chart (or *X*¯ control chart). Like the univariate counterparts, also the multivariate control charts can be distinguished into parametric and non-parametric types according to the distributive assumption underlying a control charts (e.g. normality) are verified or not. When the assumption of normality is not verified, the use of conventional (multivariate) control charts for process monitoring is questionable. Non-parametric control charts do not require distributive assumptions on process data and generally enjoy greater robustness, namely they are less sensitive to outliers than parametric control schemes. A survey of parametric multivariate SPC charts can be found in Bersimis *et al.*, 2007, while a review of non-parametric multivariate control charts can be found in Chakraborti &

Statistical depth functions are largely used in non-parametric statistics for the analysis of multivariate data. These are non-parametric functions that can provide a dimension reduction to high-dimensional problems. In this work, we will focus on *Lp* data depth to build a control chart. The *Lp* data depth have additional advantages over other existing depth based control charts already introduced in a multivariate SPC (i.e. they ensure an ease of computation even

For a more detailed description please refer to Montgomery, 2007.

ries of traditional univariate control charts (Crosier, 1988).

Graham, 2019.

in high dimensions).

Carmela Iorio 1, Giuseppe Pandolfo1, Michele Staiano1 , Massimo Aria2 and Roberta Siciliano1

1 Department of Industrial Engineering, University of Naples Federico II, Italy, (e-mail: carmela.iorio@unina, <EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
2 Department of Economics and Statistics, University of Naples Federico II, Italy<EMAIL_ADDRESS>
ABSTRACT: Control charts are used to identify non-random behaviours of a manufacturing process by monitoring changes in the distribution of the quality characteristics of the tested product. Process monitoring of related variables is usually referred to as a multivariate quality control problem. In many applications there is not enough information to justify the assumption of a specific form for the underlying process distribution. Thus, a non-parametric approach is a valid tool in a quality control process. Among possible non-parametric statistical techniques, data depth functions are gaining increasing interest in multivariate quality control. The aim of this work is to investigate the behaviour of a non-parametric approach based on the notion of the *Lp* depth in the statistical process control.

KEYWORDS: Non-parametric statistics, Q-charts, Data depth.

#### 1 Introduction

Nowadays, industries collect a large amount of data on more than one variable. Hence in a quality control process there is more than one quality variable to be monitored simultaneously. A traditional control chart monitoring a single variable is not useful for detecting the overall quality of a process, as it is determined by the interaction of several related variables (Liu, 1995, Idris *et al.*, 2019). For this reason, multivariate analysis is becoming increasingly important within the statistical process control approaches (Woodall & Montgomery, 1999). Multivariate control charts are needed when dealing with more than one quality variable as overcome the drawback of obtaining incorrect control limits when dealing with related variables. As a matter of fact, the multivariate procedure takes into account the association between the components of a multivariate process.

Multivariate quality control studies was first conducted by Hotteling, 1947. For a more detailed description please refer to Montgomery, 2007.

THE *Lp* DATA DEPTH AND ITS APPLICATION TO MULTIVARIATE PROCESS CONTROL CHARTS Carmela Iorio 1, Giuseppe Pandolfo1, Michele Staiano1 , Massimo Aria2 and Roberta Siciliano1

1 Department of Industrial Engineering, University of Naples Federico II, Italy, (e-mail: carmela.iorio@unina, <EMAIL_ADDRESS>
2 Department of Economics and Statistics, University of Naples Federico II, Italy,

ABSTRACT: Control charts are used to identify non-random behaviours of a manufacturing process by monitoring changes in the distribution of the quality characteristics of the tested product. Process monitoring of related variables is usually referred to as a multivariate quality control problem. In many applications there is not enough information to justify the assumption of a specific form for the underlying process distribution. Thus, a non-parametric approach is a valid tool in a quality control process. Among possible non-parametric statistical techniques, data depth functions are gaining increasing interest in multivariate quality control. The aim of this work is to investigate the behaviour of a non-parametric approach based on the notion of the *Lp*

Nowadays, industries collect a large amount of data on more than one variable. Hence in a quality control process there is more than one quality variable to be monitored simultaneously. A traditional control chart monitoring a single variable is not useful for detecting the overall quality of a process, as it is determined by the interaction of several related variables (Liu, 1995, Idris *et al.*, 2019). For this reason, multivariate analysis is becoming increasingly important within the statistical process control approaches (Woodall & Montgomery, 1999). Multivariate control charts are needed when dealing with more than one quality variable as overcome the drawback of obtaining incorrect control limits when dealing with related variables. As a matter of fact, the multivariate procedure takes into account the association between the components of a mul-

<EMAIL_ADDRESS><EMAIL_ADDRESS>
KEYWORDS: Non-parametric statistics, Q-charts, Data depth.

<EMAIL_ADDRESS>
depth in the statistical process control.

1 Introduction

tivariate process.

Woodall, 2000 distinguishes the techniques of the control chart processing in Phase I (also called retrospective or preliminary phase) and in Phase II. Phase I uses charts with the purpose of defining whether a process is statistically under control when the first group are processed. In Phase II, the charts are used to check if the process is in control when future subgroups were being processed. In this last phase, it is assumed that the distribution of the process is known and most of the classical applications require the hypothesis that the process under consideration follows a multivariate normal distribution. However, in most industrial applications the distribution of a parametric multivariate control chart is difficult to estimate for processes with multiple quality characteristics. As such, all observations are considered as *d*-dimensional vector and therefore used to detect possible shifts in the *d*-dimensional distributions of the quality process. A statistical process control (SPC) procedure set up in a multivariate framework is more effective than a joint monitoring system consisting of a series of traditional univariate control charts (Crosier, 1988).

The most popular multivariate statistical process control charts are based on the Hotelling's *T*2 statistics, that are a multivariate extension of Shewart's chart (or *X*¯ control chart). Like the univariate counterparts, also the multivariate control charts can be distinguished into parametric and non-parametric types according to the distributive assumption underlying a control charts (e.g. normality) are verified or not. When the assumption of normality is not verified, the use of conventional (multivariate) control charts for process monitoring is questionable. Non-parametric control charts do not require distributive assumptions on process data and generally enjoy greater robustness, namely they are less sensitive to outliers than parametric control schemes. A survey of parametric multivariate SPC charts can be found in Bersimis *et al.*, 2007, while a review of non-parametric multivariate control charts can be found in Chakraborti & Graham, 2019.

Statistical depth functions are largely used in non-parametric statistics for the analysis of multivariate data. These are non-parametric functions that can provide a dimension reduction to high-dimensional problems. In this work, we will focus on *Lp* data depth to build a control chart. The *Lp* data depth have additional advantages over other existing depth based control charts already introduced in a multivariate SPC (i.e. they ensure an ease of computation even in high dimensions).

#### 2 Our proposal

Depth-based methodology to construct control charts can be interpreted as a multivariate generalization of standard univariate. A depth function aims at providing the degree of centrality of a point *x* with respect to a distribution *F* in R*d*, denoted by *D*(*x*,*F*). Hence, higher values of *D*(*x*,*F*) correspond to deeper (more central) while smaller values indicate less central points (i.e. further away from the center with respect to *F*). Hence, a center-outward ranking of the data is provided. There are several notions of data depth function available in the literature. The halfspace, simplicial, Mahalanobis and *Lp* depths are some of the most popular ones. In this work, we adopt the notion of *Lp* depth introduced by Zuo, 2004 because of its ease of computation and (local and global) robustness properties. The depth is defined as follows:

3 Conclusion

References

In a statistical process control framework, we proposed to use control charts based on the *Lp* data depth. Our approach is fully non-parametric, meaning that the obtained charts are valid without parametric assumptions on the process distribution. In addition, these charts allow for the simultaneous detection of both the location change and the scale increase in a process. The performance of our proposal is investigated via a simulation study. The results show that the *Lp* depth based control charts are a promising alternative to the well-known Mahalanobis depth. Moreover, *Lp* depth is particularly appealing because of

BERSIMIS, SOTIRIS, PSARAKIS, STELIOS,&PANARETOS, JOHN. 2007. Multivariate statistical process control charts: an overview. *Quality and*

CHAKRABORTI, S, & GRAHAM, MA. 2019. Nonparametric (distributionfree) control charts: An updated overview and some results. *Quality En-*

CROSIER, RONALD B. 1988. Multivariate generalizations of cumulative sum

HOTTELING, H. 1947. Multivariate quality control, illustrated by the air testing of sample bombsights. *Techniques of statistical analysis*, 111–184. IDRIS, SUWANDA, WACHIDAH, LISNUR, SOFIYAYANTI, TETI,&HARA-HAP, ERWIN. 2019. The Control Chart of Data Depth Based on Influence Function of Variance Vector. *In: Journal of Physics: Conference Series*,

LIU, REGINA Y. 1995. Control charts for multivariate processes. *Journal of*

MONTGOMERY, DOUGLAS C. 2007. *Introduction to statistical quality con-*

WOODALL, WILLIAM H. 2000. Controversies and contradictions in statistical process control. *Journal of Quality Technology*, 32(4), 341–350. WOODALL, WILLIAM H, & MONTGOMERY, DOUGLAS C. 1999. Research issues and ideas in statistical process control. *Journal of Quality Technol-*

ZUO, YIJUN. 2004. Robustness of weighted L p–depth and L p–median. *All-*

*gemeines Statistisches Archiv*, 88(2), 215–234.

*the American Statistical Association*, 90(432), 1380–1387.

its computational ease even in high multidimensional spaces.

*Reliability engineering international*, 23(5), 517–543.

quality-control schemes. *Technometrics*, 30(3), 291–303.

*gineering*, 31(4), 523–544.

vol. 1366. IOP Publishing.

*trol*. John Wiley & Sons.

*ogy*, 31(4), 376–386.

$$L^p D\left(\mathbf{x}, F\right) = \frac{1}{1 + E\left(||\mathbf{x} - X||\_p\right)},$$

where *X* ∼ *F*, · denotes the *Lp*-norm (when *p* = 2 the Euclidean norm is derived) and *E* (·) is its expected value.

We conducted a Montecarlo simulation study to evaluate the performance of the *Q*-type control charts based on the *Lp* data depth in comparison with the Mahalanobis depth-based *Q*-type charts. The *Q*-type control chart is the multivariate analogue of the average univariate chart *X*¯. We set *p* = 2 for the computation of the *Lp* depth. The simulation study was designed as an analysis to evaluate the chart performances under multiple settings, defined with regard to the number of variables to be monitored (i.e., the dimension), the size of the reference sample, the size of the sub-group, and by considering different distributional settings (Normal, Skew-Normal and Cauchy). We considered both the in-control and out-of-control cases. Specifically, three out-of-control scenarios were evaluated including shift in the mean vector, change in the variance and a combination of both variance change and shift in the mean. We evaluated the performances in terms of average run length (ARL) and its standard deviation. ARL is defined as the expected number of samples required to get a first out-of-control signal, and it can be obtained by taking reciprocal of false alarm probability. Moreover, ARL is one of the performance measures used for comparing the control charts. Results obtained from both in-control and out-of-control cases indicate that *Q*-type charts based on *L*2*D* perform better than those based on Mahalanobis regardless of the process distribution, the dimensionality and the size of both the reference and the sub-group samples.

### 3 Conclusion

2 Our proposal

Depth-based methodology to construct control charts can be interpreted as a multivariate generalization of standard univariate. A depth function aims at providing the degree of centrality of a point *x* with respect to a distribution *F* in R*d*, denoted by *D*(*x*,*F*). Hence, higher values of *D*(*x*,*F*) correspond to deeper (more central) while smaller values indicate less central points (i.e. further away from the center with respect to *F*). Hence, a center-outward ranking of the data is provided. There are several notions of data depth function available in the literature. The halfspace, simplicial, Mahalanobis and *Lp* depths are some of the most popular ones. In this work, we adopt the notion of *Lp* depth introduced by Zuo, 2004 because of its ease of computation and (local and

global) robustness properties. The depth is defined as follows:

*D*(*x*,*F*) = 1

where *X* ∼ *F*, · denotes the *Lp*-norm (when *p* = 2 the Euclidean norm is

We conducted a Montecarlo simulation study to evaluate the performance of the *Q*-type control charts based on the *Lp* data depth in comparison with the Mahalanobis depth-based *Q*-type charts. The *Q*-type control chart is the multivariate analogue of the average univariate chart *X*¯. We set *p* = 2 for the computation of the *Lp* depth. The simulation study was designed as an analysis to evaluate the chart performances under multiple settings, defined with regard to the number of variables to be monitored (i.e., the dimension), the size of the reference sample, the size of the sub-group, and by considering different distributional settings (Normal, Skew-Normal and Cauchy). We considered both the in-control and out-of-control cases. Specifically, three out-of-control scenarios were evaluated including shift in the mean vector, change in the variance and a combination of both variance change and shift in the mean. We evaluated the performances in terms of average run length (ARL) and its standard deviation. ARL is defined as the expected number of samples required to get a first out-of-control signal, and it can be obtained by taking reciprocal of false alarm probability. Moreover, ARL is one of the performance measures used for comparing the control charts. Results obtained from both in-control and out-of-control cases indicate that *Q*-type charts based on *L*2*D* perform better than those based on Mahalanobis regardless of the process distribution, the dimensionality and the size of both the reference and the sub-group samples.

1+*E*(*x*−*Xp*)

,

*Lp*

derived) and *E* (·) is its expected value.

In a statistical process control framework, we proposed to use control charts based on the *Lp* data depth. Our approach is fully non-parametric, meaning that the obtained charts are valid without parametric assumptions on the process distribution. In addition, these charts allow for the simultaneous detection of both the location change and the scale increase in a process. The performance of our proposal is investigated via a simulation study. The results show that the *Lp* depth based control charts are a promising alternative to the well-known Mahalanobis depth. Moreover, *Lp* depth is particularly appealing because of its computational ease even in high multidimensional spaces.

### References

### **ANGULAR HALFSPACE DEPTH: CENTRAL REGIONS**\*

Petra Laketa1and Stanislav Nagy1

1 Faculty of Mathematics and Physics, Charles University, Prague, Czech Republic (e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
**ABSTRACT**: The angular halfspace depth is an extension of the classical halfspace depth that is applicable to directional data. The upper level sets of this function serve as an analogue of the inter-quantile regions, and allow introduction of orders and rank statistics also for data living on the unit sphere. We explore the basic theoretical properties of these regions, and contrast them with the central regions defined in multivariate Euclidean spaces using the standard halfspace depth.

**KEYWORDS**: angular depth, central regions, directional data analysis, halfspace depth.

#### **1 Angular halfspace depth**

Statistical depth is a remarkable tool for ordering multivariate data. For a probability measure *P* in the Euclidean space R*d*, *d* ≥ 1, the depth describes how much "centrally located" a point *x* ∈ R*d* is with respect to *P*. The arguably most important depth in R*d* is the *halfspace depth* that to each *x* ∈ R*d* assigns

$$hD\left(\mathbf{x}; P\right) = \inf\left\{ P\left(H\right) : H \in \mathcal{H} \text{ and } \mathbf{x} \in H \right\} \in [0, 1],\tag{1}$$

for *H* = *Hy*,*v* : *y* ∈ R*d*, *v* ∈ R*d* \ {0} the set of all closed halfspaces *Hy*,*v* = *z* ∈ R*d* : �*z*−*y*, *v*� ≥ 0 in R*d*. Here we deal with directional data (Ley & Verdebout, 2017), meaning data generated from *P* whose support lies on the unit sphere S*d*−1 = *x* ∈ R*d* : �*x*� = 1 . For most such *P*, the depth (1) is trivially zero on S*d*−1, and is therefore of no use. In that situation, it is more natural to consider an angular variant of the halfspace depth introduced by Small, 1987. Let *H*0 ⊂ *H* be the collection of those halfspaces *H*0,*v* ∈ *H* whose boundary contains the origin 0 ∈ R*d*. The *angular halfspace depth* of *x* ∈ S*d*−1 with respect to a probability measure *P* on S*d*−1 is defined as

$$ahD\left(\mathbf{x}; P\right) = \inf\left\{P\left(H\right) : H \in \mathcal{H}\_0 \text{ and } \mathbf{x} \in H\right\} \in [0, 1]. \tag{2}$$

\*This work was supported by the grant 19-16097Y of the Czech Science Foundation, and by the PRIMUS/17/SCI/3 project of Charles University. P. Laketa was supported by the OP RDE project "International mobility of research, technical and administrative staff at the Charles University" CZ.02.2.69/0.0/0.0/18 053/0016976.

The similarity of the depths (1) and (2) is obvious — one restricts to halfspaces from *H*0 when considering the angular depth. Many properties of the angular halfspace depth were explored by Liu & Singh, 1992. Here we first revise some of those known results, and then derive an array of new properties of the upper level sets of the function (2) for general probability measures on S*d*−1.

#### **2 A hemisphere of constant depth**

A peculiar property of the angular halfspace depth is that it has to be constant on a hemisphere of S*d*−1, for any *P* on S*d*−1. More precisely, Proposition 4.6 of Liu & Singh, 1992 says that there exists a hemisphere with a constant angular depth equal to α0 = inf*x*∈S*d*−1 *ahD*(*x*;*P*). That result is given without a proof and from the context it appears to be claimed for a closed hemisphere. In our first example we demonstrate that for general measures *P* one has to be cautious when formulating this statement.

Consider the probability measure *P* on the circle S1 in R2 (left panel of Figure 1) defined as a mixture of a uniform distribution on the upper halfcircle S1 + = (*x*1, *x*2) ∈ S1 : *x*2 > 0 and an atom at *a* = (1,0) with equal weights 1/2. For each *n* = 1,2,... consider a halfspace (beige region in Figure 1)

**Figure 1.** *Left: For P a mixture of uniform distribution on* S1 + *(grey arc) and an atom at point a (black point) a closed hemisphere of constant depth does not exist. Right: A measure P with five atoms such that ahD*(·;*P*) *is constant on* S1 *and equal to* α0 = 2/5*.*

$$H\_n = \left\{ (\mathbf{x}\_1, \mathbf{x}\_2) \in \mathbb{R}^2 \colon \mathbf{x}\_1 \cos(\pi/2 - 1/n) + \mathbf{x}\_2 \sin(\pi/2 - 1/n) \le \mathbf{0} \right\} \in \mathcal{H}\_0$$

not containing *a* at an angle θ*n* = −1/*n* with the *x*1-axis. Since lim*n*→∞ θ*n* = 0, surely lim*n*→∞ *P*(*Hn*) = 0, meaning that *ahD*(*x*;*P*) = 0 for any point *x* in the lower halfcircle S1 − = (*x*1, *x*2) ∈ S1 : *x*2 < 0 . We obtain α0 = 0. On the other hand, *P*(*H*) ≥ 1/2 for every *H* ∈*H*0 that contains *a*, implying that *ahD*(*a*;*P*) ≥ 1/2 > α0. Also, any *H* ∈*H*0 that contains points from S1 + is of positive *P*-mass. Overall we obtain that the angular depth is positive exactly in the set S1 + ∪ {*a*}, and there is no closed halfcircle of S1 of depth α0 = 0.

In our example there exists an open hemisphere S1 − with constant depth α0. That is not a coincidence — an analogous result for an open hemisphere of constant depth is possible to be proved\* for any measure *P* on S*d*−1. It is however interesting to note that it may also happen that the depth (2) is constant on the whole sphere S*d*−1, and therefore equal to α0 everywhere, see Figure 1.

#### **3 Central regions of the angular halfspace depth**

While for any *P* on S*d*−1 there always exists an open hemisphere *Smin* ⊂ S*d*−1 of minimal depth, its complement S*d*−1 \ *Smin* typically contains points of higher depth (2). The upper level sets of the angular halfspace depth therefore form a basis for generalizations of quantiles and inter-quantile regions to S*d*−1, in the same way as the level sets of the halfspace depth (1) do in R*d*. The *central region* of *P* at level α ≥ 0 is given by

$$D\_{\mathbf{a}} = \left\{ \mathbf{x} \in \mathbb{S}^{d-1} \colon ahD\left(\mathbf{x}; P\right) \ge \mathbf{a} \right\}.\tag{3}$$

The smallest non-empty region *D*α presents a natural analogue of the median applicable to directional data. In analogy with the corresponding properties well established for the standard halfspace depth (1) in R*d*, it is possible to show that also regions (3) posses several attractive traits — they are closed and spherically convex sets in S*d*−1 that can be represented as intersections of closed spherical halfspaces (sets of the form *H*0 ∩S*d*−1). Formal proofs of the following statements will appear in our comprehensive treatment of the theory of the angular halfspace depth that is currently in preparation. In each of the statements *P* is a Borel probability measure on S*d*−1.

**Upper semi-continuity.** The mapping S*d*−1 → [0,1]: *x* �→ *ahD*(*x*;*P*) is upper semi-continuous, i.e. for any *x* ∈ S*d*−1 and a sequence {*xn*}∞ *n*=1 ⊂ S*d*−1 that converges to *x* it holds lim sup*n*→∞ *ahD*(*xn*;*P*) ≥ *ahD*(*x*;*P*). As a consequence, all depth regions (3) are closed sets.

\*The proof of this claim is not difficult, but due to the limited available space we will present it elsewhere, together with all the other technical derivations outlined in the rest of this note.

**Intersection of halfspaces.** For any α > α0 we can write

$$D\_{\alpha} = \bigcap \left\{ \text{int}(H) : H \in \mathcal{H}\_{0} \text{ and } P(\text{int}(H)) > 1 - \alpha \right\} \cap \left( \mathbb{S}^{d-1} \backslash S\_{\text{min}} \right),$$

where int(*H*) denotes the interior of *H*, which is an open halfspace. This result is weaker than the one for the usual halfspace depth (Proposition 6 of Rousseeuw & Ruts, 1999), where one can write a central region as an intersection of closed halfspaces from *H* . As a consequence of our result we obtain that each *D*α with α > α0 is an intersection of a convex set in R*d* and a hemisphere. The case α = α0 gives trivially the whole sphere *D*α0 = S*d*−1.

#### **4 Refinements for smooth measures**

We say that a probability measure *P* on S*d*−1 is smooth if *P*(∂*H*) = 0 for any boundary hyperplane ∂*H* of *H* ∈ *H*0. It is satisfied for any *P* that has a density with respect to the spherical Lebesgue measure on S*d*−1. For smooth *P* one obtains stronger results about the central regions of *ahD*(*x*;*P*):

A useful application of these results is the construction of bagdistances for directional data presented by *H. Demni* in this book of short papers. In that contribution, bagdistances are used with success in a comprehensive simulation study of nonparametric classification of points in S*d*−1.

#### **References**

### CLUSTERING PRODUCTION INDEXES FOR CONSTRUCTION WITH FORECAST DISTRIBUTIONS

This approach accounts for the future dynamic behaviour of the time series,

*T*+*h*|*T* (·), *i* = 1,...,*S* is the forecast distribution function at a given fu-

*T*+*h*|*T* (*y*)−*F j*

*k*x*t*−*h* +*wk*0

*T*+*h*|*T* (*y*)

, conditioned on the information set available

 *r*

+*c*0 (1)

*t*−*h*),*t* = *p*+*h*,...,*T*}

*t*=*p*+*h* (*Yt* − *fmh*(x*t*−*h*;θ))2 .

*t*=*p*+*h* εˆ*t*.

*t*−*h*−*p*+1), *t* = *p* + *h*,...,*T*}, as an

*T*+*h* where ε∗

*T*+*h* is a

*T*+*h* con-

*T*−*p*−*h*+1 ∑*T*

x∗ *t*−*h*;θ 2 .

*T*+*h*|*T* is given by the law of *Y*ˆ ∗

*h* + ε∗ *dy r* = 1,2,

 *Fi*

puted directly, La Rocca *et al.*, 2021 have proposed a strategy in which the unknown distributions are consistently estimated by using a feed forward neural network estimator and the pair bootstrap approach. In particular, given the forecast horizon *h*, the unknown function *g*(·) can be approximated by using

> *ck*ψ w′

with θ = (*c*0,*c*1,..., *cm*,w1,...,w*m*,*w*10,...,*wm*0), where *m* is the hidden layer size, w*k* are the vectors of weights for the connections between input layer and hidden layer, *ck*, *k* = 1,...,*m* are the weights of the link between the hidden layer and the output; *wk*0 and *c*0 are the bias terms; ψ(·) is a proper chosen

*t*−*h* = (*Yt*−*h*,...,*Yt*−*h*−*p*+1). The general procedure is summarized in the following Algorithm.

2: Fix the hidden layer size *m* and the lag structure *p* and estimate the weights

3: Compute the residuals from the estimated network defined as: εˆ*t* = *Yt* −

*t*−*h*,...,*Y*∗

*t*=*p*+*h Y*∗ *t* − *fmh*

*YT* ,*YT*−1,...,*YT*−*p*+1;θˆ ∗

Note that, as usual, the bootstrap distribution can be approximated by Monte Carlo simulations repeating *B* times steps 5-7, and then computing the

1 *T*−*p*−*h*+1 ∑*T*

*m* ∑ *k*=1

*fmh* (x*t*−*h*;θ) =

1: Fix the forecast horizon *h* ≥ 1. Let *X* = {(*Yt*,x′

4: Compute the centered residuals: ε˜*t* = εˆ*t* − 1

*t*−*h*)=(*Y*∗

1 *T*−*p*−*h*+1 ∑*T*

random sample from the centered residuals {ε˜*t*}.

*t* ,*Y*∗

6: Get the bootstrap estimate of the neural network weights:

of the network as θˆ *h* = argminθ

*t* ,x′∗

*h* = argminθ

iid sample from the set of tuples *X*.

*T*+*h* = *fmh*

8: The bootstrap forecast distribution *F*∗

by using the *Lr*

ture point *T* +*h*, of the series y(*i*)

up to time *T*. Since the *Lr*

activation function and x′

x*t*−*h*;θˆ *h*

5: Resample {(*Y*∗

θˆ ∗

ditioned on *X*.

7: Compute *Y*ˆ ∗

where *Fi*

the network

Algorithm

*fmh*

Michele La Rocca 1, Francesco Giordano1 and Cira Perna1

1 Department of Economics and Statistics, University of Salerno, (e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
ABSTRACT: In this paper we focus on a recent proposal for clustering nonlinear time series data in which dissimilarities are computed according to time series forecast distributions. The aim is to evaluate the impact of COVID-19 pandemic on the construction sector for a set of 21 European countries.

KEYWORDS: Feedforward neural networks, bootstrap, nonlinear time series.

#### 1 Introduction

In the last decades there has been a growing interest in time series clustering. Some recent approaches rely on the use of distance criteria which compare the forecast densities estimated by using a resampling method combined with a nonparametric kernel estimator (see Alonso *et al.*, 2006 and Vilar *et al.*, 2010). More recently, La Rocca *et al.*, 2021, have proposed a novel approach for clustering nonlinear autoregressive time series based on the use of a class of neural network models to approximate the original nonlinear process, combined with the pair bootstrap as a resampling device. The aim of this paper is to discuss the novel approach and to evaluate the impact of COVID-19 pandemic on the production index for construction, an important business cycle indicator, for a set of 21 European countries.

#### 2 The clustering procedure in a nutshell

Let {*Yt*,*t* ∈ Z} be a real valued stationary stochastic process modeled as a nonlinear autoregressive (NAR) model of the form *Yt* = *g*(x*t*−1)+ε*t*, where *g*(·) is an unknown (possibly) nonlinear regression function, x′ *t*−1 = (*Yt*−1,...,*Yt*−*p*) and {ε*t*} are *iid* error terms, with E[ε*t*] = 0 and E[ε2 *t* ] > 0. Let y(1) ,...,y(*S*) be *S* observed time series of length *T* generated from a DGP of the previous class, where y(*i*) = *Y*(*i*) 1 ,...,*Y*(*i*) *T* . The aim is to cluster time series based on their full forecast distribution at a specific future time *T* +*h*, with *h* ≥ 1.

This approach accounts for the future dynamic behaviour of the time series, by using the *Lr* -norm distance *Dr*,*i j* = *Fi T*+*h*|*T* (*y*)−*F j T*+*h*|*T* (*y*) *r dy r* = 1,2, where *Fi T*+*h*|*T* (·), *i* = 1,...,*S* is the forecast distribution function at a given future point *T* +*h*, of the series y(*i*) , conditioned on the information set available up to time *T*. Since the *Lr* -norm distance previously defined cannot be computed directly, La Rocca *et al.*, 2021 have proposed a strategy in which the unknown distributions are consistently estimated by using a feed forward neural network estimator and the pair bootstrap approach. In particular, given the forecast horizon *h*, the unknown function *g*(·) can be approximated by using the network

$$f\_{mh} \left( \mathbf{x}\_{t-h}; \Theta \right) = \sum\_{k=1}^{m} c\_k \Psi \left( \mathbf{w}\_k' \mathbf{x}\_{t-h} + \mathbf{w}\_{k0} \right) + c\_0 \tag{1}$$

with θ = (*c*0,*c*1,..., *cm*,w1,...,w*m*,*w*10,...,*wm*0), where *m* is the hidden layer size, w*k* are the vectors of weights for the connections between input layer and hidden layer, *ck*, *k* = 1,...,*m* are the weights of the link between the hidden layer and the output; *wk*0 and *c*0 are the bias terms; ψ(·) is a proper chosen activation function and x′ *t*−*h* = (*Yt*−*h*,...,*Yt*−*h*−*p*+1).

The general procedure is summarized in the following Algorithm.

#### Algorithm

CLUSTERING PRODUCTION INDEXES FOR CONSTRUCTION WITH FORECAST DISTRIBUTIONS Michele La Rocca 1, Francesco Giordano1 and Cira Perna1

1 Department of Economics and Statistics, University of Salerno, (e-mail:

ABSTRACT: In this paper we focus on a recent proposal for clustering nonlinear time series data in which dissimilarities are computed according to time series forecast distributions. The aim is to evaluate the impact of COVID-19 pandemic on the con-

In the last decades there has been a growing interest in time series clustering. Some recent approaches rely on the use of distance criteria which compare the forecast densities estimated by using a resampling method combined with a nonparametric kernel estimator (see Alonso *et al.*, 2006 and Vilar *et al.*, 2010). More recently, La Rocca *et al.*, 2021, have proposed a novel approach for clustering nonlinear autoregressive time series based on the use of a class of neural network models to approximate the original nonlinear process, combined with the pair bootstrap as a resampling device. The aim of this paper is to discuss the novel approach and to evaluate the impact of COVID-19 pandemic on the production index for construction, an important business cycle indicator, for a

Let {*Yt*,*t* ∈ Z} be a real valued stationary stochastic process modeled as a nonlinear autoregressive (NAR) model of the form *Yt* = *g*(x*t*−1)+ε*t*, where *g*(·) is

be *S* observed time series of length *T* generated from a DGP of the previous

their full forecast distribution at a specific future time *T* +*h*, with *h* ≥ 1.

*t*−1 = (*Yt*−1,...,*Yt*−*p*)

,...,y(*S*)

y(1)

*t* ] > 0. Let

. The aim is to cluster time series based on

KEYWORDS: Feedforward neural networks, bootstrap, nonlinear time series.

<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
struction sector for a set of 21 European countries.

1 Introduction

set of 21 European countries.

class, where y(*i*) =

2 The clustering procedure in a nutshell

an unknown (possibly) nonlinear regression function, x′

and {ε*t*} are *iid* error terms, with E[ε*t*] = 0 and E[ε2

1 ,...,*Y*(*i*) *T*

 *Y*(*i*)

$$\hat{\boldsymbol{\Theta}}\_{h}^{\*} = \underset{\bullet}{\arg\min} \boldsymbol{\upmu}\_{\boldsymbol{\Theta}} \frac{1}{T - p\_{\boldsymbol{\bot}}^{-}h + 1} \sum\_{t=p+h}^{T} \left( Y\_{t}^{\*} - f\_{mh} \left( \mathbf{x}\_{t-h}^{\*}; \boldsymbol{\Theta} \right) \right)^{2} \dots$$

Note that, as usual, the bootstrap distribution can be approximated by Monte Carlo simulations repeating *B* times steps 5-7, and then computing the empirical cumulative distribution function (ECDF) of *Y*ˆ *b T*+*h*, *b* = 1,2,...,*B*. As a resampling device, the pair bootstrap has been implemented, a suitable choice in the context of neural network models. Moreover, being the data generating process nonlinear, a direct multi-step forecasting approach is considered, where a separate neural network model is estimated for each forecasting horizon, and forecasts are computed only conditioning on the observed data.

#### 3 An application to the European construction sector

The proposed procedure has been used to cluster the production index for construction (seasonally and calendar adjusted) for 21 European countries observed from January 2000 to December 2020 (base year 2015). The production index measures the activity in the building and construction industry, and it is considered a critical business cycle indicator. The dataset is available from the Eurostat website. The aim here is to identify the different group structure induced by the COVID-19 pandemic by using the forecast one-step ahead distribution for January 2020 (so excluding any observations from the COVID-19 pandemic), the forecast twelve-step ahead distribution for January 2021, the forecast one-step ahead distribution for January 2021 (we have trained all models up to December 2020).

Hungary

Hungary Bulgaria Slovakia Slovenia Croatia Netherlands Austria Sweden Poland Romania Finland Luxembourg Germany Denmark Czechia UK Spain Italy France Belgium Portugal

Cluster Dendrogram

UK Czechia Portugal Bulgaria Slovakia Spain Belgium France Hungary Denmark Germany Slovenia Croatia Netherlands Romania Poland Austria Sweden Italy Luxembourg Finland

tributions and *L*1-norm distance.

0

10

Height

20

0

10

20

Height

30

40

0

10

20

Height

30

40

Bulgaria

Slovakia

Cluster Dendrogram

Cluster Dendrogram

Spain

Italy

Portugal

Belgium

France

Romania

UK

(a) Training period January 2000 - December 2019, prediction *h* = 1

(b) Training period January 2000 - December 2019, prediction *h* = 12

(c) Training period January 2000 - December 2020, prediction *h* = 1 Figure 1: Construction indexes clustering based on *h*-step ahead forecast dis-

Czechia

Luxembourg

Germany

Finland

Denmark

Poland

Slovenia

Sweden

Netherlands

Croatia

Austria

Apparently, the group structure would have been almost identical without the impact due to the COVID-19 pandemic, showing a somewhat stable economic evolution of all the countries considered in the application (see panels a and b). On the contrary, when based on models that include the year 2020 in the training period, where all countries experienced severe contractions in their economic activities, the dataset shows a pretty different group structure, indicating different routes and timelines for economic recovery (see panel c).

#### References

empirical cumulative distribution function (ECDF) of *Y*ˆ *b*

3 An application to the European construction sector

models up to December 2020).

References

*ted*.

As a resampling device, the pair bootstrap has been implemented, a suitable choice in the context of neural network models. Moreover, being the data generating process nonlinear, a direct multi-step forecasting approach is considered, where a separate neural network model is estimated for each forecasting horizon, and forecasts are computed only conditioning on the observed data.

The proposed procedure has been used to cluster the production index for construction (seasonally and calendar adjusted) for 21 European countries observed from January 2000 to December 2020 (base year 2015). The production index measures the activity in the building and construction industry, and it is considered a critical business cycle indicator. The dataset is available from the Eurostat website. The aim here is to identify the different group structure induced by the COVID-19 pandemic by using the forecast one-step ahead distribution for January 2020 (so excluding any observations from the COVID-19 pandemic), the forecast twelve-step ahead distribution for January 2021, the forecast one-step ahead distribution for January 2021 (we have trained all

Apparently, the group structure would have been almost identical without the impact due to the COVID-19 pandemic, showing a somewhat stable economic evolution of all the countries considered in the application (see panels a and b). On the contrary, when based on models that include the year 2020 in the training period, where all countries experienced severe contractions in their economic activities, the dataset shows a pretty different group structure, indicating different routes and timelines for economic recovery (see panel c).

ALONSO, A.M., BERRENDERO, J.R., HERNANDEZ ´ , A., & JUSTEL, A. 2006. Time series clustering based on forecast densities. *Computational*

LA ROCCA, M., GIORDANO, F., & PERNA, C. 2021. Clustering nonlinear time series with neural network bootstrap forecast distributions. *Submit-*

VILAR, J.A., ALONSO, A.M., & VILAR, J.M. 2010. Non-linear time series clustering based on non-parametric forecast densities. *Computational*

*Statistics & Data Analysis*, 51(2), 762–776.

*Statistics & Data Analysis*, 54(11), 2850–2865.

*T*+*h*, *b* = 1,2,...,*B*.

(a) Training period January 2000 - December 2019, prediction *h* = 1 Cluster Dendrogram

(b) Training period January 2000 - December 2019, prediction *h* = 12 Cluster Dendrogram

(c) Training period January 2000 - December 2020, prediction *h* = 1 Figure 1: Construction indexes clustering based on *h*-step ahead forecast distributions and *L*1-norm distance.

## CLUSTERING LONGITUDINAL DATA WITH CATEGORY THEORY FOR DIABETIC KIDNEY DISEASE

ory\* (Grandis, 2020), we can exploit its concepts for cluster analysis (Carlsson & Mémoli, 2013). We introduce comparisons between patient trajectories according to their shapes. Here, we compare patients' trajectories building clusters of shapes using the Fréchet distance (Genolini *et al.*, 2016), which takes into account shape variations. The novelty of our study is an integrative approach joining categories, clusters, and shape trajectories. For each patient we observe demographic and clinical variables, treatments, and response to the different treatments. This approach is derived for a real dataset concerning patients with diabetic kidney disease (DKD) from the DC-ren project.† Clustering is based on the response to the treatment, that is evaluated from the estimated glomerular filtration rate (eGFR). The identification of different evolution patterns can shed light on the best individualized drug combination. The paper is structured as follows: in Section 2 we introduce some theoretical concepts and the methodology we adopted, and in Section 3 we analyze the

Let us consider a dataset composed of *n* patients characterized by *p* observable variables at four time points *t*0, *t*1, *t*2, *t*3. Each patient is characterized as a triplet (X*i*(*tk*), D(*tk*), *Yi*(*tk*)), where *i* is the individual (the patient); *tk* is the time point *k* = 0,1,2,3; X*i*(*tk*) is a set of variables characterizing the individual; *Yi*(*tk*) is the value of the response variable *Y* at *tk*; D(*tk*) stands for the

obtain an *enriched double category* with metrics in R (Grandis, 2020), having

The comparison of trajectories of different patients involves both of these dis-

of significative curves, which groups similar trajectories in the same cluster.

\*A category is constituted by objects (points) and morphisms (arrows). A *functor* maps objects and morphisms of a category into objects and morphisms of another category. *Natural*

with respect to the

*i* (*tk*,*tk*). We thus

*i* . Clustering

*i* (*tk*,*tk*) as morphisms.

(*tk*), and the distance between values observed at

(*tk*), *dY*

ι >, ι = 1,...,*n* < *n*, where *n* is the number

*i*,*i*

given drug. We indicate the distance between patients *i*,*i*

times *tk*,*tk* of the variable *Y* for the same individual *i* as *dY*

tances. The time trajectory of the *i*-th patient is indicated as *pathY*

†https://dc-ren.eu/. The project focuses on type 2 diabetes.

*i*,*i*

*i* →< *patha*,*Y*

*i* (*tk*), *k* = 0,...3, *i* = 1,...,*n*, as objects, and *dY*

results of our study.

2 Methodology

variable *Y* and time *k* as *dY*

is a functor *Fa* : *pathY*

*transformations* map functors to functors.

*xY*

Maria Mannone12, Veronica Distefano1, Claudio Silvestri13 and Irene Poli1

1 European Centre for Living Technology, Ca' Foscari University of Venice, Italy, (e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
2 Department of Mathematics and Computer Sciences, University of Palermo, Italy

3 Dipartimento di Scienze Ambientali, Informatica e Statistica, Ca' Foscari University of Venice, Italy

ABSTRACT: In the framework of precision medicine, we investigate the similarity of diabetic kidney disease (DKD) patients through longitudinal data clusters. Starting with insights from category theory, we build patients' clusters according to the shapes of their trajectories, adopting the Fréchet distance. We group patients according to their behavior of the estimated glomerular filtration rate (eGFR), obtaining informative mean curves. Behavior pattern recognition can shed light on individualized treatments.

KEYWORDS: longitudinal data clustering, category theory, Fréchet distance, precision medicine, DKD disease progress

### 1 Introduction

Precision medicine aims to find individualized therapeutic treatments according to patients' specific characteristics. To make accurate predictions it is crucial to retrieve information on the long-term reactions of patients to given treatments, investigating time trajectories of the disease progress (Karpati *et al.*, 2018). Each patient is identified by demographic and clinical variables at different time points. The similarity of behavior of patients across time can be accounted by *clusters of trajectories*. The final aim of this research is to build clusters of longitudinal data to identify the optimal treatment rule. Therefore, we intend to build patient clusters according to distances between patients at each time point, and distances of the same patient between time points. We can consider distance as a transformation, and distance variation as a transformation between transformations. Because the concept of transformations between transformations is the starting point of mathematical category theory\* (Grandis, 2020), we can exploit its concepts for cluster analysis (Carlsson & Mémoli, 2013). We introduce comparisons between patient trajectories according to their shapes. Here, we compare patients' trajectories building clusters of shapes using the Fréchet distance (Genolini *et al.*, 2016), which takes into account shape variations. The novelty of our study is an integrative approach joining categories, clusters, and shape trajectories. For each patient we observe demographic and clinical variables, treatments, and response to the different treatments. This approach is derived for a real dataset concerning patients with diabetic kidney disease (DKD) from the DC-ren project.† Clustering is based on the response to the treatment, that is evaluated from the estimated glomerular filtration rate (eGFR). The identification of different evolution patterns can shed light on the best individualized drug combination. The paper is structured as follows: in Section 2 we introduce some theoretical concepts and the methodology we adopted, and in Section 3 we analyze the results of our study.

#### 2 Methodology

CLUSTERING LONGITUDINAL DATA WITH CATEGORY THEORY FOR DIABETIC KIDNEY DISEASE Maria Mannone12, Veronica Distefano1, Claudio Silvestri13 and Irene Poli1

1 European Centre for Living Technology, Ca' Foscari University of Venice, Italy, (e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
2 Department of Mathematics and Computer Sciences, University of Palermo, Italy 3 Dipartimento di Scienze Ambientali, Informatica e Statistica, Ca' Foscari University

ABSTRACT: In the framework of precision medicine, we investigate the similarity of diabetic kidney disease (DKD) patients through longitudinal data clusters. Starting with insights from category theory, we build patients' clusters according to the shapes of their trajectories, adopting the Fréchet distance. We group patients according to their behavior of the estimated glomerular filtration rate (eGFR), obtaining informative mean curves. Behavior pattern recognition can shed light on individualized treatments.

KEYWORDS: longitudinal data clustering, category theory, Fréchet distance, preci-

Precision medicine aims to find individualized therapeutic treatments according to patients' specific characteristics. To make accurate predictions it is crucial to retrieve information on the long-term reactions of patients to given treatments, investigating time trajectories of the disease progress (Karpati *et al.*, 2018). Each patient is identified by demographic and clinical variables at different time points. The similarity of behavior of patients across time can be accounted by *clusters of trajectories*. The final aim of this research is to build clusters of longitudinal data to identify the optimal treatment rule. Therefore, we intend to build patient clusters according to distances between patients at each time point, and distances of the same patient between time points. We can consider distance as a transformation, and distance variation as a transformation between transformations. Because the concept of transformations between transformations is the starting point of mathematical category the-

<EMAIL_ADDRESS><EMAIL_ADDRESS>
of Venice, Italy

1 Introduction

sion medicine, DKD disease progress

Let us consider a dataset composed of *n* patients characterized by *p* observable variables at four time points *t*0, *t*1, *t*2, *t*3. Each patient is characterized as a triplet (X*i*(*tk*), D(*tk*), *Yi*(*tk*)), where *i* is the individual (the patient); *tk* is the time point *k* = 0,1,2,3; X*i*(*tk*) is a set of variables characterizing the individual; *Yi*(*tk*) is the value of the response variable *Y* at *tk*; D(*tk*) stands for the given drug. We indicate the distance between patients *i*,*i* with respect to the variable *Y* and time *k* as *dY i*,*i* (*tk*), and the distance between values observed at times *tk*,*tk* of the variable *Y* for the same individual *i* as *dY i* (*tk*,*tk*). We thus obtain an *enriched double category* with metrics in R (Grandis, 2020), having *xY i* (*tk*), *k* = 0,...3, *i* = 1,...,*n*, as objects, and *dY i*,*i* (*tk*), *dY i* (*tk*,*tk*) as morphisms. The comparison of trajectories of different patients involves both of these distances. The time trajectory of the *i*-th patient is indicated as *pathY i* . Clustering is a functor *Fa* : *pathY i* →< *patha*,*Y* ι >, ι = 1,...,*n* < *n*, where *n* is the number of significative curves, which groups similar trajectories in the same cluster.

\*A category is constituted by objects (points) and morphisms (arrows). A *functor* maps objects and morphisms of a category into objects and morphisms of another category. *Natural transformations* map functors to functors.

†https://dc-ren.eu/. The project focuses on type 2 diabetes.

The < *patha*,*Y* ι >, ∀ι is the representative curve of each cluster.‡ Most of the existing research uses the Euclidean distance to compare trajectories. However, because we aim to compare trajectory shapes, we determine the Fréchet distance, and, according to this distance, we build patient clusters. The Fréchet distance is based on the comparison between pairs of points following the profiles of the curves they belong to. We analyzed the response to the treatment, measured according to the eGFR variable. We build patient clusters deriving the mean of eGFR trajectories. We then investigate the characteristics of patients, in relationship with demographic and clinical variables which characterize each patient. We evaluated the behavior of 241 DKD patients observed in a 4-year period, according to the DC-ren project. Computationally, we used an extension of the longitudinal k-means, *kmlShape* (Genolini *et al.*, 2016), with time scale 0.5.

do not have a positive response. In cluster 7, all patients are keeping the same drug (*D*1), and most of them show a positive response. Patients in clusters 6 and 7 display the best eGFR behavior; most of the patients in these clusters keep in fact a stable behavior and a positive response to the treatment. Patients in clusters 1 and 6 start from close eGFR values, but different UACR values; given *D*1 +*D*2, their response is different. Patients in clusters 1, 2, 5 change, receiving *D*1 +*D*3 and *D*1 +*D*4. These results will be considered in building a

Figure 1: Shape clusters (a) and Table of mean values (b).

This research activity is part of the project DC-ren that has received funding from the European Union's Horizon 2020 research and innovation program under grant agreement No. 848011. We thank the researchers of the European Centre for Living Technology (ECLT) for helpful discussions and suggestions.

CARLSSON, G., & MÉMOLI, F. 2013. Classifying Clustering Schemes. *Foun-*

GENOLINI, C., ECOCHARD, R., BENGHEZAL, M., DRISS, T., ANDRIEU, S., & SUBTIL, F. 2016. kmlShape: An Efficient Method to Cluster Longitudinal Data (Time-Series) According to Their Shapes. *PlosOne*, 11(6). GRANDIS, M. 2020. *Higher Category Theory*. Singapore: World Scientific. KARPATI, T., LEVENTER-ROBERTS, M., FELDMAN, B., C., COHEN-STAVI., & RAZ I., BALICER, R. 2018. Patient clusters based on HbA1c trajectories: A step toward individualized medicine in type 2 diabetes.

*dations of Computational Mathematics*, 13, 221–252.

**variables / clusters cl. 1 (19%) cl. 2 (16%) cl. 3 (13%) cl 4 (12%) cl. 5 (11%) cl. 6 (11%) cl. 7 (10%) cl. 8 (7%)** eGFR (t0) 82 (6) 71 (5) 67 (7) 44 (8) 49 (6) 79 (9) 57 (6) 37 (4) eGFR (t1) 79 (11) 68 (9) 57 (11) 43 (7) 52 (6) 89 (14) 62 (9) 31 (6) eGFR (t2) 76 (10) 67 (9) 52 (10) 39 (6) 51 (5) 83 (10) 64 (9) 29 (7) eGFR (t3) 69 (7) 63 (6) 46 (7) 39 (6) 51 (5) 89 (10) 66 (6) 26 (4) mean\_UACR (t0) 29.36 (48.07) 58.27 (278.21) 170.18 (466.31) 162.48 (513.48) 46.76 (143.67) 40.41 (85.99) 46.41 (144.17) 142.71 (246.52) mean\_UACR (t1) 22.61 (25.89) 51.99 (162.43) 84.76 (296.70) 85.46 (208.09) 42.50 (136.90) 57.50 (200.63) 36.38 (108.26) 126.52 (220.99) mean\_UACR (t2) 25.15 (13.29) 56.10 (219.01) 100.46 (254.98) 112.26 (236.83) 82.67 (227.14) 31.24 (71.27) 39.98 (81.98) 131.71 (300.80) mean\_UACR (t3) 47.07 (120.65) 40.55 (130.90) 121.71 (337.43) 107.78 (262.89) 76.96 (177.44) 29.66 (65.78) 71.67 (205.72) 219.31 (553.02) HbA1c (t0) 7.0 (1.3) 7.2 (1.3) 7.2 (1.1) 7.1 (1.3) 7.2 (1.1) 7.4 (0.9) 7.1 (1.3) 6.8 (0.8) HbA1c (t1) 7.2 (1.2) 7.5 (1.3) 7.3 (1.2) 7.2 (1.7) 7.3 (1.2) 7.3 (1.2) 7.5 (1.3) 7.0 (1.4) HbA1c (t2) 7.1 (1.2) 7.4 (1.5) 7.2 (1.1) 7.4 (1.3) 7.2 (1.0) 7.2 (1.0) 7.3 (1.2) 7.3 (1.7) HbA1c (t3) 7.0 (1.0) 7.4 (1.2) 7.4 (1.1) 7.8 (1.4) 7.9 (1.7) 7.9 (1.7) 7.5 (1.2) 7.9 (1.9) (b)

predictive system to envisage the best treatment for each individual.

0.5 **cl. 1**

1.0 1.5 2.0 2.5 3.0 3.5 4.0

t0 t1 t2 t3 time (a)

7% **cl. 2 cl. 3 cl. 4 cl. 5 cl. 6 cl. 7 cl. 8**

19% 16% 13% 12% 11% 11% 10%

**cl. 1 cl. 6 cl. 2 cl. 3 cl. 7 cl. 5 cl. 4 cl. 8**

Times 1.0 1.5 2.0 2.5 3.0 3.5 4.0

Time

Times

*PLos One*, 13(11).

Acknowledgements

References

eGFR

#### 3 Results

In this study, which involves clusters based on shapes of individual trajectories, without assuming a particular shape, we obtained a grouping of patients according to their similarity of eGFR behavior. Trajectories are evaluated upon the Fréchet distance between them, computed on the continuous variable eGFR. We obtain 8 patient clusters with similar eGFR shape of individual trajectories. In Figure 1a, we represent the obtained clusters. In the Table (Figure 1b), we show the behavior of clusters achieved with the Fréchet distance, with standard deviations at each time point. We find three main patients' subgroups: patients with decreasing eGFR (clusters 1, 3, 8); low decreasing eGFR (cl. 2, 4), and stable/increasing eGFR (cl. 5, 6, 7). To explain these behaviors, we consider some relevant clinical characteristics of patients, which are the ratio of urinary albumin and creatinine (mean UACR) and the glycated hemoglobin (HbA1c), shown in the Table (Figure 1b). Mean UACR is decreasing in cluster 6, while it is increasing in cluster 8. Decreasing or stable values of HbA1c, which characterize DKD patients (Karpati *et al.*, 2018), are shown by patients in cl. 6, while increasing values of HbA1c are shown by patients in cl. 8. We notice that there is a relationship between non-decreasing eGFR and stable HbA1c. Patients can receive different treatments, such as *D*1, *D*2, *D*3, *D*4. In cluster 3 most of the patients that change drug (*D*1 → *D*1 +*D*2) show a positive response to the treatment. However, patients in cluster 2 who change the drug

‡A different clustering method *Fb* gives us similar representative paths < *pathb*,*Y* ι >. Natural transformation α*a*,*b* : *Fa* → *Fb* maps clustering methods.

do not have a positive response. In cluster 7, all patients are keeping the same drug (*D*1), and most of them show a positive response. Patients in clusters 6 and 7 display the best eGFR behavior; most of the patients in these clusters keep in fact a stable behavior and a positive response to the treatment. Patients in clusters 1 and 6 start from close eGFR values, but different UACR values; given *D*1 +*D*2, their response is different. Patients in clusters 1, 2, 5 change, receiving *D*1 +*D*3 and *D*1 +*D*4. These results will be considered in building a predictive system to envisage the best treatment for each individual.

Times Figure 1: Shape clusters (a) and Table of mean values (b).

#### Acknowledgements

This research activity is part of the project DC-ren that has received funding from the European Union's Horizon 2020 research and innovation program under grant agreement No. 848011. We thank the researchers of the European Centre for Living Technology (ECLT) for helpful discussions and suggestions.

#### References

The < *patha*,*Y*

with time scale 0.5.

3 Results

ι >, ∀ι is the representative curve of each cluster.‡ Most of the

existing research uses the Euclidean distance to compare trajectories. However, because we aim to compare trajectory shapes, we determine the Fréchet distance, and, according to this distance, we build patient clusters. The Fréchet distance is based on the comparison between pairs of points following the profiles of the curves they belong to. We analyzed the response to the treatment, measured according to the eGFR variable. We build patient clusters deriving the mean of eGFR trajectories. We then investigate the characteristics of patients, in relationship with demographic and clinical variables which characterize each patient. We evaluated the behavior of 241 DKD patients observed in a 4-year period, according to the DC-ren project. Computationally, we used an extension of the longitudinal k-means, *kmlShape* (Genolini *et al.*, 2016),

In this study, which involves clusters based on shapes of individual trajectories, without assuming a particular shape, we obtained a grouping of patients according to their similarity of eGFR behavior. Trajectories are evaluated upon the Fréchet distance between them, computed on the continuous variable eGFR. We obtain 8 patient clusters with similar eGFR shape of individual trajectories. In Figure 1a, we represent the obtained clusters. In the Table (Figure 1b), we show the behavior of clusters achieved with the Fréchet distance, with standard deviations at each time point. We find three main patients' subgroups: patients with decreasing eGFR (clusters 1, 3, 8); low decreasing eGFR (cl. 2, 4), and stable/increasing eGFR (cl. 5, 6, 7). To explain these behaviors, we consider some relevant clinical characteristics of patients, which are the ratio of urinary albumin and creatinine (mean UACR) and the glycated hemoglobin (HbA1c), shown in the Table (Figure 1b). Mean UACR is decreasing in cluster 6, while it is increasing in cluster 8. Decreasing or stable values of HbA1c, which characterize DKD patients (Karpati *et al.*, 2018), are shown by patients in cl. 6, while increasing values of HbA1c are shown by patients in cl. 8. We notice that there is a relationship between non-decreasing eGFR and stable HbA1c. Patients can receive different treatments, such as *D*1, *D*2, *D*3, *D*4. In cluster 3 most of the patients that change drug (*D*1 → *D*1 +*D*2) show a positive response to the treatment. However, patients in cluster 2 who change the drug

‡A different clustering method *Fb* gives us similar representative paths < *pathb*,*Y*

transformation α*a*,*b* : *Fa* → *Fb* maps clustering methods.

ι >. Natural

GRANDIS, M. 2020. *Higher Category Theory*. Singapore: World Scientific. KARPATI, T., LEVENTER-ROBERTS, M., FELDMAN, B., C., COHEN-

STAVI., & RAZ I., BALICER, R. 2018. Patient clusters based on HbA1c trajectories: A step toward individualized medicine in type 2 diabetes. *PLos One*, 13(11).

## A REDUNDANCY ANALYSIS WITH MULTIVARIATE RANDOM-COEFFICIENTS LINEAR MODELS

*Analysis* (*RDA)* was originally introduced in order to capture the effect onto a reduced space of the linear dependence by a set of criterion variables from a set of predictors. A *RDA* of the predicted criterion variables by the best linear unbiased predictor may be quite representative (Marcis & Salvatore, 2020). This paper uses a *RDA* by a least-squares solution for an optimal fixed-effects estimate from the data provided by the random-coefficients linear model predictors of the criterion variables. The application study performs the method introduced on the official data by the Italian Equitable and Sustainable Well-

2 Redundancy Analysis: model estimation and application study

Given a *q*-variate random vector *Y*, consider the case when *Y* is partitioned in *n* subjects (groups), each of them with *ni* individuals (*i* = 1,...,*m*,; *j* = 1,...,*ni*;*s* = 1,...,*q*). We assume that the population model for the *n* sub-

gression coefficients. *Ai* is a matrix of *q*-variate *r*-dimensional vectors of

*r* = *p*, the population model is a multivariate *RCM*, with *zi* = *xi*. Given a sample of *N* = Σ*i*Σ*jni j* units (e.g. repeated measurements), then the model struc-

covariates, *Z*+ the design matrix of random effects, and *E* the matrix of regression within-subject errors, *cov*(*vec*(*E*)) = *R*. Assuming both *Y* and *X*+ as columnwise centered and standardized, we get in the general RCM setup

*y*∗, β = *vec*(*B*), Σ = *var*(*y*) = *E*

tains the row joint reduced coordinates in the space of *Y*∗∗. In accordance with recent law reforms, the Equitable and Sustainable Well-being indicators (BES) - annually provided by the Italian Statistical Institute(ISTAT, 2017) are designed to define the economic policies which largely act on some fundamental aspects of the quality of life. In order to highlight the result of the proposed *RDA,* we use 12 BES indicators relating to the years 2013-2016, collected at NUTS-2 level. In particular the variables are S8 (Age-standardised mortality rate), IF3 (People with tertiary education), L12 (Satisfaction with

=*Y* = *X*+

*Y* , further noticing that *U*∗

*si*. If *F* = *Y*∗∗*var*(*y*)−1

*Y* ∗∗

*i*

*q*×*pxi*|*q*×1 + *Ai*|*q*×*rzi*|*r*×1, where *B* is the matrix of fixed re-

) ∼ *N*(0,Σ*a*), Σ*a* = *cov*(*vec*(*A*

*i*,*ss*)), where *s*,*s* = 1,...,*q*, the *r*×*r* blocks of Σ*a* . When

*N*×*pmApm*×*q* + *EN*×*q*, with *X*+ the matrix of data

 (*y*∗ *s* −*y*∗ *s*) (*y*∗ *s* −*y*∗ *s*) ,

*i*

2 − *X*+*B*,*Y*∗∗ = *Y* − *E*(*Y*), and

β gives the common rescaled

*Y*∗∗ = *YV* −1Λ*Y* con-

)) = {Σ*a*,*ss*},

being indicators.

jects is *yi*|*q*×1 = *B*

ture is *YN*×*q* = *X*+

*cov*(*asi*, *ysi*) = *DZ*

Σ−1

*X*)−1*X*

predictor's coordinates, *UY*Λ*YV*

 β = (*X*

and Σ*a*,*ss* = *cov*(*vec*(*A*

random-effects, with *ai* = *vec*(*A*

*N*×*pBp*×*q* + *Z*+

*i* = *DX*

Σ−1

the singular value decomposition of

Laura Marcis 1, Maria Chiara Pagliarella 2 and Renato Salvatore 1

1 Department of Economics and Law, University of Cassino, (e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>2 Italian National Institute for Public Policy Analysis (e-mail<EMAIL_ADDRESS>
ABSTRACT: Random-coefficients linear models can be considered as a particular case of linear mixed models, in which the random effects depend on the model fixedeffects design matrix. A Redundancy Analysis of estimates of the multivariate randomeffects may be able to capture the leading contribution to this correlation. Starting from the standardized multivariate best linear predictors, we introduce the random effects reduced space by a weighted least-squares closed-form solution. The application shows the effect of the linear dependence of the random-effects in the space of the predictor variables.

KEYWORDS: Redundancy analysis, linear mixed model, empirical best linear unbiased predictor, restricted maximum likelihood estimator.

### 1 Introduction

Random-coefficients linear regression models (*RCM*) represent a special case of linear mixed models (LMM, Demidenko, 2004), where the vector of regression coefficients for the subjects (e.g. repeated observations) is modeled in a second stage linear regression equation. In order to specify this type of models, it is convenient to define a two-stage hierarchical linear model, with a first stage that models within-subject observations, and as second stage we use a linear model for the random regression coefficients. Although in the basic linear mixed models the random effects are not correlated with the modeled response variables (unlike the fixed-effects with the random effects estimates), in the *RCM* this correlation depends on the fixed-effects design matrix of the regression model. Since this happens, one can be interested to know in which components of a multivariate model the random effects are related to the subspace spanned by the model covariates. In the same way, what components of the multivariate vector seem to be orthogonal to that subspace. *Redundancy* *Analysis* (*RDA)* was originally introduced in order to capture the effect onto a reduced space of the linear dependence by a set of criterion variables from a set of predictors. A *RDA* of the predicted criterion variables by the best linear unbiased predictor may be quite representative (Marcis & Salvatore, 2020). This paper uses a *RDA* by a least-squares solution for an optimal fixed-effects estimate from the data provided by the random-coefficients linear model predictors of the criterion variables. The application study performs the method introduced on the official data by the Italian Equitable and Sustainable Wellbeing indicators.

A REDUNDANCY ANALYSIS WITH MULTIVARIATE RANDOM-COEFFICIENTS LINEAR MODELS Laura Marcis 1, Maria Chiara Pagliarella 2 and Renato Salvatore 1

1 Department of Economics and Law, University of Cassino, (e-mail:

2 Italian National Institute for Public Policy Analysis (e-mail:

ABSTRACT: Random-coefficients linear models can be considered as a particular case of linear mixed models, in which the random effects depend on the model fixedeffects design matrix. A Redundancy Analysis of estimates of the multivariate randomeffects may be able to capture the leading contribution to this correlation. Starting from the standardized multivariate best linear predictors, we introduce the random effects reduced space by a weighted least-squares closed-form solution. The application shows the effect of the linear dependence of the random-effects in the space of the

KEYWORDS: Redundancy analysis, linear mixed model, empirical best linear unbi-

Random-coefficients linear regression models (*RCM*) represent a special case of linear mixed models (LMM, Demidenko, 2004), where the vector of regression coefficients for the subjects (e.g. repeated observations) is modeled in a second stage linear regression equation. In order to specify this type of models, it is convenient to define a two-stage hierarchical linear model, with a first stage that models within-subject observations, and as second stage we use a linear model for the random regression coefficients. Although in the basic linear mixed models the random effects are not correlated with the modeled response variables (unlike the fixed-effects with the random effects estimates), in the *RCM* this correlation depends on the fixed-effects design matrix of the regression model. Since this happens, one can be interested to know in which components of a multivariate model the random effects are related to the subspace spanned by the model covariates. In the same way, what components of the multivariate vector seem to be orthogonal to that subspace. *Redundancy*

<EMAIL_ADDRESS><EMAIL_ADDRESS>
ased predictor, restricted maximum likelihood estimator.

<EMAIL_ADDRESS>
predictor variables.

1 Introduction

#### 2 Redundancy Analysis: model estimation and application study

Given a *q*-variate random vector *Y*, consider the case when *Y* is partitioned in *n* subjects (groups), each of them with *ni* individuals (*i* = 1,...,*m*,; *j* = 1,...,*ni*;*s* = 1,...,*q*). We assume that the population model for the *n* subjects is *yi*|*q*×1 = *B q*×*pxi*|*q*×1 + *Ai*|*q*×*rzi*|*r*×1, where *B* is the matrix of fixed regression coefficients. *Ai* is a matrix of *q*-variate *r*-dimensional vectors of random-effects, with *ai* = *vec*(*A i* ) ∼ *N*(0,Σ*a*), Σ*a* = *cov*(*vec*(*A i* )) = {Σ*a*,*ss*}, and Σ*a*,*ss* = *cov*(*vec*(*A i*,*ss*)), where *s*,*s* = 1,...,*q*, the *r*×*r* blocks of Σ*a* . When *r* = *p*, the population model is a multivariate *RCM*, with *zi* = *xi*. Given a sample of *N* = Σ*i*Σ*jni j* units (e.g. repeated measurements), then the model structure is *YN*×*q* = *X*+ *N*×*pBp*×*q* + *Z*+ *N*×*pmApm*×*q* + *EN*×*q*, with *X*+ the matrix of data

covariates, *Z*+ the design matrix of random effects, and *E* the matrix of regression within-subject errors, *cov*(*vec*(*E*)) = *R*. Assuming both *Y* and *X*+ as columnwise centered and standardized, we get in the general RCM setup *cov*(*asi*, *ysi*) = *DZ i* = *DX si*. If *F* = *Y*∗∗*var*(*y*)−1 2 − *X*+*B*,*Y*∗∗ = *Y* − *E*(*Y*), and β = (*X* Σ−1 *X*)−1*X* Σ−1 *y*∗, β = *vec*(*B*), Σ = *var*(*y*) = *E* (*y*∗ *s* −*y*∗ *s*) (*y*∗ *s* −*y*∗ *s*) , the singular value decomposition of *Y* ∗∗ =*Y* = *X*+ β gives the common rescaled predictor's coordinates, *UY*Λ*YV Y* , further noticing that *U*∗ *Y*∗∗ = *YV* −1Λ*Y* contains the row joint reduced coordinates in the space of *Y*∗∗. In accordance with recent law reforms, the Equitable and Sustainable Well-being indicators (BES) - annually provided by the Italian Statistical Institute(ISTAT, 2017) are designed to define the economic policies which largely act on some fundamental aspects of the quality of life. In order to highlight the result of the proposed *RDA,* we use 12 BES indicators relating to the years 2013-2016, collected at NUTS-2 level. In particular the variables are S8 (Age-standardised mortality rate), IF3 (People with tertiary education), L12 (Satisfaction with

job), REL4 (Social participation), POL5 (Trust in institutions), SIC1 (Homicide rate), BS3 (Positive judgment for future perspectives), PATR9 (Presence of Parks/Gardens), AMB9 (Satisfaction for life), INN1 (Percentage of R&D expenditure), Q2 (Childhood services) and LBE1 (logarithm of per-capita adjusted disposable income). We use the latter as the predictor variable in the *RCM*, while the remaining 11 variables are dependent variables. The application uses the restricted maximum likelihood estimation, inside a SAS/IML code. To simplify the estimation process, we assume equicorrelation between the multivariate components of random effects. The linear mixed model with random coefficients highlights its analytical capabilities in the Figure 1. The plot features the standardized best predictors (STDP) and the original criterion variables in the space of the latter. As an example, while there is no correlation between INN1 (R&D expenditure) and AMB9 (satisfaction for the environment), the same best predictor variables register an inverse correlation. This evidence is supported by arguments, such as critical differences among Northern and Southern Italian Regions. Figure 2 shows the constrained RDA, in which the major contribution is given by the variables IF3, Q2, REL4, and AMB9. Interpreting the correlation between random effects and the LBE1 model covariate, this dependence is mainly explained indeed by the amount of the population that completed tertiary education (IF3). This means that most of the differences between Italian Regions reflect the dependence of the IF3 variable on the disposable income.

References

Sons.

DEMIDENKO, E. 2004. *Mixed Models: theory and applications*. Wiley and

MARCIS, L., & SALVATORE, R. 2020. Joint Redundancy Analysis by a Multivariate Linear Predictor. *Conference of the Italian Statistical Society*.

#### References

job), REL4 (Social participation), POL5 (Trust in institutions), SIC1 (Homicide rate), BS3 (Positive judgment for future perspectives), PATR9 (Presence of Parks/Gardens), AMB9 (Satisfaction for life), INN1 (Percentage of R&D expenditure), Q2 (Childhood services) and LBE1 (logarithm of per-capita adjusted disposable income). We use the latter as the predictor variable in the *RCM*, while the remaining 11 variables are dependent variables. The application uses the restricted maximum likelihood estimation, inside a SAS/IML code. To simplify the estimation process, we assume equicorrelation between the multivariate components of random effects. The linear mixed model with random coefficients highlights its analytical capabilities in the Figure 1. The plot features the standardized best predictors (STDP) and the original criterion variables in the space of the latter. As an example, while there is no correlation between INN1 (R&D expenditure) and AMB9 (satisfaction for the environment), the same best predictor variables register an inverse correlation. This evidence is supported by arguments, such as critical differences among Northern and Southern Italian Regions. Figure 2 shows the constrained RDA, in which the major contribution is given by the variables IF3, Q2, REL4, and AMB9. Interpreting the correlation between random effects and the LBE1 model covariate, this dependence is mainly explained indeed by the amount of the population that completed tertiary education (IF3). This means that most of the differences between Italian Regions reflect the dependence of the IF3

variable on the disposable income.

### THE USE OF MULTIPLE IMPUTATION TECHNIQUES IN SOCIAL MEDIA DATA

step, the substitution of missing values is implemented using a threshold based on the number of expressed "Likes". In particular, a missing value is considered as a "Dislike" only when a user has expressed a percentage of "Likes" that is higher than a selected threshold. Alternatively, if the percent of "Likes" is less than the threshold, a missing "Like" is imputed as a "Nothing". The second step has been pursued using a multiple imputation technique known as MIMCA method (Multiple Imputation with Multiple Correspondence Analysis) (Audigier *et al.*, 2017). This procedure is applied to social media data from

Multiple Imputation with Multiple Correspondence Analysis represents an available alternative as imputation technique for qualitative data. This approach allows to impute data sets with incomplete categorical variables. The principle of MI with MCA, as well as all the other multiple imputation techniques, consists in creating *M* different datasets to reflect the uncertainty on imputed values. In this context, each dataset is obtained with an algorithm called *iterative MCA*, which is useful to impute qualitative data. The iterative MCA algorithm consists in recoding the incomplete dataset as an incomplete disjunctive table *Z*, randomly imputing the missing values, estimating the principal components and loadings from the completed matrix and then, using these estimates to impute missing values according to the following reconstruction formula:

*Z*ˆ = *U*ˆΛˆ*V*ˆ *T* + *M*.

where *U*ˆ, Λˆ and *V*ˆ are the left singular vectors, the diagonal matrix of singular values and the right singular vectors, respectively. The final version of matrix *Z* is obtained as *Z* = *W* ∗ *Z* + (I −*W*) ∗ *Z*ˆ, where ∗ is the Hadamard product and *W* is a matrix of weights where *wi j* = 1 if *zi j* is missing and *wi j* = 0 otherwise. In this context, MCA is configured as a singular values decomposition

the disjunctive table, *M* is a matrix whose rows are equal to the vectors of the means of each component of *Z*, *D*Σ is a diagonal matrix with the proportions of individuals characterized by a specific category and *R* is the matrix of uniform weights assigned to individuals. After a first step of imputation, the procedure of iterative MCA is repeated many times until a convergence criterion is reached. In many cases, due to overfitting problems, a regularized version of

*K D*−1

Σ ,*R*). The matrix *Z* represents

the official pages of 7 Italian newspapers.

applied on the triplet of matrices (*Z* − *M*, 1

this algorithm is used (Josse *et al.*, 2012).

2 The MIMCA approach

Paolo Mariani 1, Andrea Marletta 1 and Matteo Locci 1

1 Department of Economics, Management and Statistics, University of Milano-Bicocca, (e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>)

ABSTRACT: In the big data context, it is very frequent to manage the analysis of missing values. This is especially relevant in the field of statistical analysis, where this represents a thorny issue. This study proposes a strategy for data enrichment in presence of sparse matrices. The research objective consists in the evaluation of a possible distinction of behaviour among observations in sparse matrices with missing data. After selecting among the multiple imputation methods, an innovative technique will be presented to impute missing observations as a negative position or a neutral opinion. This method has been applied to a dataset measuring the interaction between users and social network pages for some Italian newspapers.

KEYWORDS: Social network data, Missing values, Multiple imputations

#### 1 Introduction

The treatment of missing values is still a neglected phase in the field of quantitative analysis. In not statistical contexts, the most abused solution is the row elimination, that is to say the deletion of the observation with missing values. This operation could result misleading and the treatment of missing observations is more complex procedure. Firstly, it is necessary to conduct some preliminary analysis about the nature of this lack of information and to recognize the mechanism of the missing data. This relationship aims to evaluate the link between the observed value and the missing one. This lead to the well-known classification of Little and Rubin (Little, 1988) in MCAR (Missing Completely At Random) data, MAR (Missing at Random) data o NMAR (Not Missing At Random) data. Only after the identification of these mechanism, it is possible to find the best solution to solve the problem of missing values. If the complete case analysis has not been considered as a valid alternative, it is necessary to proceed with the imputation of the missing observation.

In this study, an innovative technique to discern the missing value from a behaviour for some individuals has been proposed using two steps. In the first step, the substitution of missing values is implemented using a threshold based on the number of expressed "Likes". In particular, a missing value is considered as a "Dislike" only when a user has expressed a percentage of "Likes" that is higher than a selected threshold. Alternatively, if the percent of "Likes" is less than the threshold, a missing "Like" is imputed as a "Nothing". The second step has been pursued using a multiple imputation technique known as MIMCA method (Multiple Imputation with Multiple Correspondence Analysis) (Audigier *et al.*, 2017). This procedure is applied to social media data from the official pages of 7 Italian newspapers.

#### 2 The MIMCA approach

THE USE OF MULTIPLE IMPUTATION TECHNIQUES IN SOCIAL MEDIA DATA Paolo Mariani 1, Andrea Marletta 1 and Matteo Locci 1

1 Department of Economics, Management and Statistics, University of Milano-Bicocca, (e-mail<EMAIL_ADDRESS>
ABSTRACT: In the big data context, it is very frequent to manage the analysis of missing values. This is especially relevant in the field of statistical analysis, where this represents a thorny issue. This study proposes a strategy for data enrichment in presence of sparse matrices. The research objective consists in the evaluation of a possible distinction of behaviour among observations in sparse matrices with missing data. After selecting among the multiple imputation methods, an innovative technique will be presented to impute missing observations as a negative position or a neutral opinion. This method has been applied to a dataset measuring the interaction between

<EMAIL_ADDRESS><EMAIL_ADDRESS>)

users and social network pages for some Italian newspapers.

proceed with the imputation of the missing observation.

1 Introduction

KEYWORDS: Social network data, Missing values, Multiple imputations

The treatment of missing values is still a neglected phase in the field of quantitative analysis. In not statistical contexts, the most abused solution is the row elimination, that is to say the deletion of the observation with missing values. This operation could result misleading and the treatment of missing observations is more complex procedure. Firstly, it is necessary to conduct some preliminary analysis about the nature of this lack of information and to recognize the mechanism of the missing data. This relationship aims to evaluate the link between the observed value and the missing one. This lead to the well-known classification of Little and Rubin (Little, 1988) in MCAR (Missing Completely At Random) data, MAR (Missing at Random) data o NMAR (Not Missing At Random) data. Only after the identification of these mechanism, it is possible to find the best solution to solve the problem of missing values. If the complete case analysis has not been considered as a valid alternative, it is necessary to

In this study, an innovative technique to discern the missing value from a behaviour for some individuals has been proposed using two steps. In the first Multiple Imputation with Multiple Correspondence Analysis represents an available alternative as imputation technique for qualitative data. This approach allows to impute data sets with incomplete categorical variables. The principle of MI with MCA, as well as all the other multiple imputation techniques, consists in creating *M* different datasets to reflect the uncertainty on imputed values. In this context, each dataset is obtained with an algorithm called *iterative MCA*, which is useful to impute qualitative data. The iterative MCA algorithm consists in recoding the incomplete dataset as an incomplete disjunctive table *Z*, randomly imputing the missing values, estimating the principal components and loadings from the completed matrix and then, using these estimates to impute missing values according to the following reconstruction formula:

$$
\hat{\mathbf{Z}} = \hat{U}\hat{\Lambda}\hat{\mathcal{V}}^T + M.
$$

where *U*ˆ, Λˆ and *V*ˆ are the left singular vectors, the diagonal matrix of singular values and the right singular vectors, respectively. The final version of matrix *Z* is obtained as *Z* = *W* ∗ *Z* + (I −*W*) ∗ *Z*ˆ, where ∗ is the Hadamard product and *W* is a matrix of weights where *wi j* = 1 if *zi j* is missing and *wi j* = 0 otherwise. In this context, MCA is configured as a singular values decomposition applied on the triplet of matrices (*Z* − *M*, 1 *K D*−1 Σ ,*R*). The matrix *Z* represents the disjunctive table, *M* is a matrix whose rows are equal to the vectors of the means of each component of *Z*, *D*Σ is a diagonal matrix with the proportions of individuals characterized by a specific category and *R* is the matrix of uniform weights assigned to individuals. After a first step of imputation, the procedure of iterative MCA is repeated many times until a convergence criterion is reached. In many cases, due to overfitting problems, a regularized version of this algorithm is used (Josse *et al.*, 2012).

This approach is part of the family of joint modelling MI method, which means that it is more computationally efficient than conditional models. In fact, this MI technique is based on Multiple Correspondence Analysis and then the number of parameters estimated is small. Another advantage of MI with MCA is the goodness of estimation even if the number of individuals is small. Finally, MI with MCA well represents less frequent categories in the step of imputation. This last is another property that derives from MCA.

the following:

References

*and computing*, 27(2), 501–518.

83(404), 1198–1202.

*Journal of classification*, 29(1), 91–116.

then a "Dislike" is imputed;

"Dislike" nor a "Nothing" is imputed.

a "Nothing" is imputed;

• if the proportion of "Dislike" imputed is greater than or equal to 60%,

• if the proportion of "Dislike" imputed is less than or equal to 40%, then

• if the proportion of "Dislike" is between 40% and 60%, then neither a

La Repubblica Corriere della Sera Il Fatto Quotidiano Il Sole 24 Ore

La Gazzetta dello Sport Il Messaggero La Stampa Total

As can be noted from the table 1, few missing values are still present. In fact, in some cases the number of "Dislike" imputed in a specific cell is very similar to the number of "Nothing". In particular, this behaviour is manifested when the proportion of "Dislike" (or "Nothing") is between 40% and 60%. Even if there are some cases where missing values are not imputed, MI with MCA works well. In fact, the proportion of missing values is now equal to 6%.

AUDIGIER, V, HUSSON, F, & JOSSE, J. 2017. MIMCA: multiple imputation for categorical variables with multiple correspondence analysis. *Statistics*

JOSSE, J., CHAVENT, M., LIQUET, B., & HUSSON, F. 2012. Handling missing values with regularized iterative multiple correspondence analysis.

LITTLE, R. 1988. A test of missing completely at random for multivariate data with missing values. *Journal of the American statistical Association*,

Table 1. *Distribution of "Like", "Dislike" and "Nothing" after MI with MCA.*

"Like" 299 244 268 158 "Dislike" 24 8 57 47 "Nothing" 149 235 139 197 Missing Values 24 9 32 94

"Like" 116 66 86 1237 "Dislike" 221 255 200 812 "Nothing" 155<PHONE_NUMBER> Missing Values 4 6 40 209

#### 3 Application on Italian newspaper social pages

The dataset used for this application is represented by users that expressed at least one "Like" in social media pages, websites, and forums concerning drugs and health. The research was conducted on 2,795 Italian subjects considering all interactions between people and brands and between products and services on Facebook. The selected category for Facebook pages is Italian newspapers. Each column of the dataset is a dummy variable that represents the presence or absence of a "Like." The 7 Italian newspapers are: La Repubblica, Corriere della Sera, Il Fatto Quotidiano, Il Sole 24 Ore, La Gazzetta dello Sport, Il Messaggero, La Stampa.

Before performing the MIMCA approach, the entire dataset has been divided into training and validation set. In particular, a number of cells equivalent to 30% of the cells observed has been set to "missing value". In order to create a validation set similar to the original data set, the proportion of each category ("Like", "Dislike" and "Nothing") has been maintained. The number of multiple data sets generated is equal to 100. The category to be imputed is selected by the majority rule. In other terms, among 100 imputations for each cell of the validation set, the category imputed at least 34 times is selected.

In order to evaluate the performances of MI with MCA, a confusion matrix has been created and summarized through an index of accuracy. This approach imputes more than 81% of the cells considered. Then, the performance of this technique is satisfactory.

Once the goodness of MIMCA has been proved, the process of data enrichment about cells without a category observed or imputed can be completed. In order to achieve this goal, *M* = 100 data sets have been imputed with MIMCA and, for each cell with a missing value, only those where a "Dislike" or a "Nothing" has been imputed are considered. Moreover, in order to minimize the simulation error due to the application of a bootstrap procedure, a threshold has been introduced. More specifically, considering only the data sets where a "Dislike" or a "Nothing" has been imputed, the imputation rule for each cell is the following:

This approach is part of the family of joint modelling MI method, which means that it is more computationally efficient than conditional models. In fact, this MI technique is based on Multiple Correspondence Analysis and then the number of parameters estimated is small. Another advantage of MI with MCA is the goodness of estimation even if the number of individuals is small. Finally, MI with MCA well represents less frequent categories in the step of

The dataset used for this application is represented by users that expressed at least one "Like" in social media pages, websites, and forums concerning drugs and health. The research was conducted on 2,795 Italian subjects considering all interactions between people and brands and between products and services on Facebook. The selected category for Facebook pages is Italian newspapers. Each column of the dataset is a dummy variable that represents the presence or absence of a "Like." The 7 Italian newspapers are: La Repubblica, Corriere della Sera, Il Fatto Quotidiano, Il Sole 24 Ore, La Gazzetta dello Sport, Il

Before performing the MIMCA approach, the entire dataset has been divided into training and validation set. In particular, a number of cells equivalent to 30% of the cells observed has been set to "missing value". In order to create a validation set similar to the original data set, the proportion of each category ("Like", "Dislike" and "Nothing") has been maintained. The number of multiple data sets generated is equal to 100. The category to be imputed is selected by the majority rule. In other terms, among 100 imputations for each cell of the validation set, the category imputed at least 34 times is selected.

In order to evaluate the performances of MI with MCA, a confusion matrix has been created and summarized through an index of accuracy. This approach imputes more than 81% of the cells considered. Then, the performance of this

Once the goodness of MIMCA has been proved, the process of data enrichment about cells without a category observed or imputed can be completed. In order to achieve this goal, *M* = 100 data sets have been imputed with MIMCA and, for each cell with a missing value, only those where a "Dislike" or a "Nothing" has been imputed are considered. Moreover, in order to minimize the simulation error due to the application of a bootstrap procedure, a threshold has been introduced. More specifically, considering only the data sets where a "Dislike" or a "Nothing" has been imputed, the imputation rule for each cell is

imputation. This last is another property that derives from MCA.

3 Application on Italian newspaper social pages

Messaggero, La Stampa.

technique is satisfactory.

Table 1. *Distribution of "Like", "Dislike" and "Nothing" after MI with MCA.*

As can be noted from the table 1, few missing values are still present. In fact, in some cases the number of "Dislike" imputed in a specific cell is very similar to the number of "Nothing". In particular, this behaviour is manifested when the proportion of "Dislike" (or "Nothing") is between 40% and 60%. Even if there are some cases where missing values are not imputed, MI with MCA works well. In fact, the proportion of missing values is now equal to 6%.

#### References

### PREDICTION OF GENE EXPRESSION FROM TRANSCRIPTION FACTORS AFFINITIES: AN APPLICATION OF BAYESIAN NON-LINEAR MODELLING

Person log Expr. of EGFR

*EE* Response (log y)

Figure 1. *Left: Humans have two copies of DNA in each cell; the expression of a gene is the amount of RNA it produces. Transcription factors bind the DNA at the regulatory regions from where they activate or*

into numbers to be used in a regression model. Most existing methods rely on genotypes (discrete variables taking values 0, 1, or 2 encoding single-letter differences in the DNA of different people), which do not allow for easy interpretation (*e.g.*, "If the DNA has an 'A' instead of a 'T', the expression of the gene will be higher"). Our first goal is to develop a more interpretable model. Gene expression is mainly controlled by specialised proteins called transcription factors, which bind the DNA at particular locations (regulatory regions) by establishing weak chemical bonds. Different DNA sequences will have, therefore, different chemical *affinities* for the transcription factors. Since different individuals have different DNA sequences, it is possible to use the affinities for transcription factors as numerical (continuous) predictors in the predictive model of gene expression. Affinities have a far superior interpretation, exemplified by statements such as "If the affinity for this transcription

However, one needs to make an assumption about the relationship (*e.g.*, linear) between affinities and gene expression. de Boer *et al.* , 2020 models the logarithm of the expression as a linear function of the affinities. The model is developed for a type of yeast and achieves a good performance, but is still too simple for our application. Indeed, yeast has two important distinguishing features: 1) it is haploid, meaning that it has only one copy of DNA, whereas humans have two; and 2) its genes are regulated primarily by one regulatory

In this paper, we set up a predictive model for the expression of the EGFR (Epidermal Growth Factor Receptor) gene, and explicitly address both limitations in de Boer *et al.* , 2020. Figure 1 provides a schematic of our application.

*inhibit the expression of their target gene. Right: A plausible instance of our data set.*

factor is higher, the expression will be higher."

region, whereas human genes typically have more than one.

*AAAAAAAAAAAAAAAAAA* Predictors (A)

Affinity of TF *l* for copy 2 of region *r*

Federico Marotta1, Paolo Provero1 and Silvia Montagna2,3

1 Dipartimento di Neuroscienze "Rita Levi Montalcini", Universita` degli Studi di Torino, Via Cherasco, 15, 10126, Torino, Italy (e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
2 Dipartimento di Scienze Economico-sociali e Matematico-statistiche, Universita` degli Studi di Torino, Corso Unione Sovietica, 218/bis, 10134 Torino, Italy, (e-mail<EMAIL_ADDRESS>
3 Collegio Carlo Alberto, Piazza Vincenzo Arbarello, 8, 10122 Torino, Italy

ABSTRACT: The prediction of gene expressions from DNA sequences is a relevant problem in biology. While most of the existing methods dedicated to this task use genotypes as predictors, here we propose a method based on transcription factor affinities, which have a clearer biological interpretation. This novelty, however, introduces new challenges for modelling, which we address leveraging on Bayesian non-linear modelling techniques.

KEYWORDS: Bayesian Methods, Gene Expressions, Non-linear Predictive Modelling.

#### 1 Introduction

Scientists are often interested in predicting differences in the expression of a gene in different individuals solely from the DNA sequence of the individuals. The predicted expression can then be used in place of the real one when measuring the latter is too expensive, and the learnt relationship between DNA and expression can lead to a better understanding of how genes are regulated (Manor & Segal, 2013). The expression of a gene is the amount of RNA molecules it produces. Humans have two independent sets of DNA molecules, one coming from the father and one from the mother, therefore there are two copies of each gene. When measuring the expression, one simply sums the molecules produced by each copy.

When associating DNA to gene expressions, the first problem we face is how to encode the DNA (a 3-billion letter string from the alphabet {*A*,*C*,*G*,*T*})

PREDICTION OF GENE EXPRESSION FROM TRANSCRIPTION FACTORS AFFINITIES: AN APPLICATION OF BAYESIAN NON-LINEAR MODELLING Federico Marotta1, Paolo Provero1 and Silvia Montagna2,3

1 Dipartimento di Neuroscienze "Rita Levi Montalcini", Universita` degli Studi di Torino, Via Cherasco, 15, 10126, Torino, Italy (e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
2 Dipartimento di Scienze Economico-sociali e Matematico-statistiche, Universita` degli Studi di Torino, Corso Unione Sovietica, 218/bis, 10134 Torino, Italy, (e-mail:

ABSTRACT: The prediction of gene expressions from DNA sequences is a relevant problem in biology. While most of the existing methods dedicated to this task use genotypes as predictors, here we propose a method based on transcription factor affinities, which have a clearer biological interpretation. This novelty, however, introduces new challenges for modelling, which we address leveraging on Bayesian non-linear

KEYWORDS: Bayesian Methods, Gene Expressions, Non-linear Predictive Modelling.

Scientists are often interested in predicting differences in the expression of a gene in different individuals solely from the DNA sequence of the individuals. The predicted expression can then be used in place of the real one when measuring the latter is too expensive, and the learnt relationship between DNA and expression can lead to a better understanding of how genes are regulated (Manor & Segal, 2013). The expression of a gene is the amount of RNA molecules it produces. Humans have two independent sets of DNA molecules, one coming from the father and one from the mother, therefore there are two copies of each gene. When measuring the expression, one simply sums the

When associating DNA to gene expressions, the first problem we face is how to encode the DNA (a 3-billion letter string from the alphabet {*A*,*C*,*G*,*T*})

3 Collegio Carlo Alberto, Piazza Vincenzo Arbarello, 8, 10122 Torino, Italy

<EMAIL_ADDRESS>
modelling techniques.

1 Introduction

molecules produced by each copy.

Figure 1. *Left: Humans have two copies of DNA in each cell; the expression of a gene is the amount of RNA it produces. Transcription factors bind the DNA at the regulatory regions from where they activate or inhibit the expression of their target gene. Right: A plausible instance of our data set.*

into numbers to be used in a regression model. Most existing methods rely on genotypes (discrete variables taking values 0, 1, or 2 encoding single-letter differences in the DNA of different people), which do not allow for easy interpretation (*e.g.*, "If the DNA has an 'A' instead of a 'T', the expression of the gene will be higher"). Our first goal is to develop a more interpretable model.

Gene expression is mainly controlled by specialised proteins called transcription factors, which bind the DNA at particular locations (regulatory regions) by establishing weak chemical bonds. Different DNA sequences will have, therefore, different chemical *affinities* for the transcription factors. Since different individuals have different DNA sequences, it is possible to use the affinities for transcription factors as numerical (continuous) predictors in the predictive model of gene expression. Affinities have a far superior interpretation, exemplified by statements such as "If the affinity for this transcription factor is higher, the expression will be higher."

However, one needs to make an assumption about the relationship (*e.g.*, linear) between affinities and gene expression. de Boer *et al.* , 2020 models the logarithm of the expression as a linear function of the affinities. The model is developed for a type of yeast and achieves a good performance, but is still too simple for our application. Indeed, yeast has two important distinguishing features: 1) it is haploid, meaning that it has only one copy of DNA, whereas humans have two; and 2) its genes are regulated primarily by one regulatory region, whereas human genes typically have more than one.

In this paper, we set up a predictive model for the expression of the EGFR (Epidermal Growth Factor Receptor) gene, and explicitly address both limitations in de Boer *et al.* , 2020. Figure 1 provides a schematic of our application.

#### 2 Methodology and results

Our dataset consists in the expression values of the EGFR gene for *n* = 414 individuals (from The GTEx Consortium, 2020), and in the affinity of each regulatory region for all transcription factors, for a total of *p* = 358 predictors. Table 1. *Results of the nested-cross validation. MSE is the mean squared error,* ρ *the correlation between true and predicted expression; averages and standard deviations of these quantities are computed across the 5-folds. Avg R*2 *is the average of the squared correlations. Z is the Z-score computed via Stouffer's*

Gene Avg MSE Sd MSE Avg ρ Sd ρ Avg *R*2 *Z* pval *Z* EGFR 0.011 0.003 0.199 0.065 0.043 4.030 2.8e-5

estimator of β when imposing a Gaussian prior on β. 2) Ability to encode knowledge via priors distributions. For example, we can exploit existing biological data about which transcription factors are bound to a region of DNA by

To carry out an unbiased evaluation of the performance, we implemented a 5-fold cross-validation strategy. Table 1 summarises the results. While the average *R*2 may seem small, we emphasise that low values are common in the prediction of gene expression and our model outperforms recently published genotype-based models (the *R*2 achieved by Nagpal *et al.* , 2019 is only 0.005). Thus, our method can model the underlying biological problem in a realistic way and provide meaningful results thanks to its interpretable predictors. In the future, it could be improved by considering interactions between transcription factors, which are also biologically important. Nevertheless, for the time being, we hope that non-linear models will find their way in the field of gene expression prediction, which currently is dominated by genotype-based

DE BOER, CARL G., VAISHNAV, EESHIT D., SADEH, RONEN, *et al.* . 2020. Deciphering eukaryotic gene-regulatory logic with 100 million random

MANOR, OHAD,&SEGAL, ERAN. 2013. Robust Prediction of Expression Differences among Human Individuals Using Only Genotype Informa-

NAGPAL, SINI, MENG, XIAORAN, EPSTEIN, MICHAEL P., *et al.* . 2019. TIGAR: An Improved Bayesian Tool for Transcriptomic Data Imputation Enhances Gene Mapping of Complex Traits. *Am. J. Hum. Genet.*, 105(2),

THE GTEX CONSORTIUM. 2020. The GTEx Consortium atlas of genetic regulatory effects across human tissues. *Science*, 369(6509), 1318–1330.

promoters. *Nat. Biotechnol.*, 38(1), 56–65.

tion. *PLoS Genet.*, 9(3), e1003396.

*method, which combines the* ρ *of the five folds, and pval Z is the p-value of the Z-score.*

giving the corresponding variables a less stringent regularisation.

linear models.

References

258–266.

We can take multiple regulatory regions into account (goal 2 above) via a straightforward modification of the model in de Boer *et al.* , 2020, which becomes: log(*y*) = β0 + ∑*r g*=1 ∑*l f*=1 A*f g*β*f g*. Here *y* denotes the gene expression, {*Af g*}*l*,*r f*=1,*g*=1 are the affinities, and β = (β0,...,β*rl*) is a vector of model parameters. Similarly to de Boer *et al.* , 2020, we sum over all transcriptions factors, indexed by *f* , but now also along the regulatory regions, *g*, of the gene.

Accommodating for both copies of DNA (goal 1 above) is more challenging. Biologically, we know that the effects of the two copies should be additive in the original scale of the expression, not in the log-transformed expression. At the same time, working with the expression in the original scale can be troublesome, for it is often not normally distributed. Therefore, we propose the following model for the expression of a single gene:

$$\log(\mathfrak{y}) \sim \text{mvnormal}\left(\log\left(e^{\mathbf{A}^{(1)}\mathfrak{B}} + e^{\mathbf{A}^{(2)}\mathfrak{B}}\right), \sigma^2 I\right). \tag{1}$$

Here *y* is an *n*-vector of expression values (one for each individual), *A*(*i*) , with *i* ∈ {1,2}, is the *n* × *rl* affinity matrix for copy *i*, where each column represents a transcription factor-regulatory region pair (*r* is the number of regions, *l* the number of transcription factors), and vector β (*lr*×1) encapsulates the coefficients of the affinities. By computing the exponential of *A*(*i*) β, with *i* ∈ {1,2}, we obtain the effect of copy *i* on the expression in the original scale. We subsequently sum the two effects, and take the log of the sum to go back to the log-scale response. Importantly, the coefficient of a given transcription factor in a given regulatory region is the same for the two copies of DNA. We notice that for this reason our model does not fall in the class of generalised linear models (at least not obviously), as each coefficient β*j* appears two times independently for two different predictors.

Model (1) is embedded in a Bayesian framework by placing a normal prior (with mean zero and variance τ) on all coefficients β independently, and a noninformative Jeffreys prior on σ2. The Bayesian framework is chosen for two primary reasons: 1) Although in our specific application *p* < *n*, generally in genomics *p n* and regularisation is needed. Regularisation can be thought of as imposing a Bayesian prior on the underlying parameters. Specifically, the ridge regression estimator can be viewed as the Bayesian posterior mean

Table 1. *Results of the nested-cross validation. MSE is the mean squared error,* ρ *the correlation between true and predicted expression; averages and standard deviations of these quantities are computed across the 5-folds. Avg R*2 *is the average of the squared correlations. Z is the Z-score computed via Stouffer's method, which combines the* ρ *of the five folds, and pval Z is the p-value of the Z-score.*

estimator of β when imposing a Gaussian prior on β. 2) Ability to encode knowledge via priors distributions. For example, we can exploit existing biological data about which transcription factors are bound to a region of DNA by giving the corresponding variables a less stringent regularisation.

To carry out an unbiased evaluation of the performance, we implemented a 5-fold cross-validation strategy. Table 1 summarises the results. While the average *R*2 may seem small, we emphasise that low values are common in the prediction of gene expression and our model outperforms recently published genotype-based models (the *R*2 achieved by Nagpal *et al.* , 2019 is only 0.005).

Thus, our method can model the underlying biological problem in a realistic way and provide meaningful results thanks to its interpretable predictors. In the future, it could be improved by considering interactions between transcription factors, which are also biologically important. Nevertheless, for the time being, we hope that non-linear models will find their way in the field of gene expression prediction, which currently is dominated by genotype-based linear models.

#### References

2 Methodology and results

becomes: log(*y*) = β0 + ∑*r*

sion, {*Af g*}*l*,*r*

Our dataset consists in the expression values of the EGFR gene for *n* = 414 individuals (from The GTEx Consortium, 2020), and in the affinity of each regulatory region for all transcription factors, for a total of *p* = 358 predictors. We can take multiple regulatory regions into account (goal 2 above) via a straightforward modification of the model in de Boer *et al.* , 2020, which

parameters. Similarly to de Boer *et al.* , 2020, we sum over all transcriptions factors, indexed by *f* , but now also along the regulatory regions, *g*, of the gene. Accommodating for both copies of DNA (goal 1 above) is more challenging. Biologically, we know that the effects of the two copies should be additive in the original scale of the expression, not in the log-transformed expression. At the same time, working with the expression in the original scale can be troublesome, for it is often not normally distributed. Therefore, we propose

*f*=1,*g*=1 are the affinities, and β = (β0,...,β*rl*) is a vector of model

β +*eA*(2) β ,σ2 *I*

. (1)

,

β, with

*f*=1 A*f g*β*f g*. Here *y* denotes the gene expres-

*g*=1 ∑*l*

the following model for the expression of a single gene:

 log *eA*(1)

Here *y* is an *n*-vector of expression values (one for each individual), *A*(*i*)

with *i* ∈ {1,2}, is the *n* × *rl* affinity matrix for copy *i*, where each column represents a transcription factor-regulatory region pair (*r* is the number of regions, *l* the number of transcription factors), and vector β (*lr*×1) encapsulates the coefficients of the affinities. By computing the exponential of *A*(*i*)

*i* ∈ {1,2}, we obtain the effect of copy *i* on the expression in the original scale. We subsequently sum the two effects, and take the log of the sum to go back to the log-scale response. Importantly, the coefficient of a given transcription factor in a given regulatory region is the same for the two copies of DNA. We notice that for this reason our model does not fall in the class of generalised linear models (at least not obviously), as each coefficient β*j* appears two times

Model (1) is embedded in a Bayesian framework by placing a normal prior (with mean zero and variance τ) on all coefficients β independently, and a noninformative Jeffreys prior on σ2. The Bayesian framework is chosen for two primary reasons: 1) Although in our specific application *p* < *n*, generally in genomics *p n* and regularisation is needed. Regularisation can be thought of as imposing a Bayesian prior on the underlying parameters. Specifically, the ridge regression estimator can be viewed as the Bayesian posterior mean

log(*y*) ∼ mvnormal

independently for two different predictors.

## HIGH DIMENSIONAL MODEL-BASED CLUSTERING OF EUROPEAN GEOREFERENCED VEGETATION PLOTS

on the Western Australian continental margin (Woolley et al., 2013) and forest physiognomic types in Italy (Attorre et al., 2014). However, when we face with high-dimensional vegetation data, FMM, or more specifically, standard model-based clustering techniques, may show a disappointing behavior. This is mainly due to the fact that the number of parameters to be estimated usually depends on the dimension of the observed space and such approaches may therefore suffer from the so-called curse of dimensionality (Bellman, 1957). In this paper, we suggest the use of a robust model-based clustering, named Gaussian mixture models for high-dimensional data (HD-GMM) proposed by Bouveyron et al. (2007). We examine a database of 7955 georeferenced plots and 3181 plant species of evergreen forest vegetation, created in TURBOVEG by storing published and unpublished phytosociological plots collected over the last 30 years. These plant communities are scattered along the Mediterranean coastal area, whose main and distinctive ecological feature is the prolonged aridity in summer and rainfall mainly concentrated during winter and spring. Making use of HD-GMM, we assume that high-dimensional vetegation data live in subspaces with a dimensionality lower than the dimensionality of the original plant species space, limiting the number of parameters to estimate and, consequently, the computational time. Finally, as in the FMM framework, the plots are classified based on their plant species composition through a posteriori specific-plot probability and the clusters are defined to be homogeneous

in that they include plots that show similar vegetation.

Let Y = (Y1,...,Y*n*) be the abundance data matrix, where the generic element *yi j* represents the value of the measure of abundance for the *j*-th tree species, namely the percentage of biomass of a certain species with respect to the total biomass of vegetation, observed in the *i*-th plot of the study area, (*i* = 1,...,*n*, *j* = 1,..., *p*). FMM assumes that each plot y*i* is drawn from a mixture of *K* components in some unknown mixing proportions π1,...,π*K*,

*k*=1 π*k* = 1. Each component identifies a cluster. When a (multivariate) Gaussian density is used to describe the component-specific distribution of observed plant species cover, the component is identified by a specific center, defined by the mean vector (as the observed values are on abundance scale, we may hypothesize that similar plots will be characterized by similar values of abundance of the same species), and a specific shape, summarized by the covariance matrix, which allows for varying dependence between cover values corresponding to different plant species for plots in that component. In other

2 The model

with ∑*K*

Francesca Martella 1, Fabio Attorre 2, Michele De Sanctis 2 and Giuliano Fanelli 2

1 Department of Statistical Sciences, Sapienza University of Rome, (e-mail<EMAIL_ADDRESS>
2 Department of Environmental Biology, Sapienza University of Rome, Italy (e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
ABSTRACT: An important challenge in complex vegetation systems is the classification of vegetation since it represents a useful tool for summarizing our knowledge of vegetation patterns and, consequently, for nature conservation, landscape mapping and land-use planning. It typically requires standard clustering methods that are capable of identifying groups of plots characterized by dominant and diagnostic species. When the data are high-dimensional, however, efficient clustering methods have to be considered. In this paper, we consider a robust model-based clustering, called Gaussian mixture models for high-dimensional data (HD-GMM) which takes into account for the specific subspace around which each cluster is located and, consequently, provides parsimonious modeling. Results are encouraging and deserve further discussion.

KEYWORDS: vegetation plots, high-dimensional data, finite mixture models

#### 1 Introduction