Cristian Preda

Quentin Grimonprez

Vincent Vandewalle

cfda: an R Package for Categorical Functional Data Analysis

Keywords: functional data, categorical data, stochastic process, multiple correspondence analysis

Categorical functional data represented by paths of a stochastic jump process with continuous time and finite set of states are considered. As an extension of the multiple correspondence analysis to an infinite set of variables, optimal encodings of states over time are approximated using an arbitrary finite basis of functions. That allows dimension reduction, optimal representation and visualisation of data in lower dimensional spaces. The methodology is implemented in the cfda R package and is illustrated using a real data set in the clustering framework.

Introduction

Most literature devoted to functional data considers data as sample paths of a real-valued stochastic process, X = {X t , t ∈ T }, X t ∈ R p , p ≥ 1 where T is some continuous set. Among a considerable record of papers on the subject, the monographs of [START_REF] Ramsay | Functional Data Analysis[END_REF] and [START_REF] Ferraty | Nonparameric Functional Data Analysis[END_REF] are still the main references presenting methodologies for visualisation, denoising, clustering and regression when dealing with functional data represented by real-valued functions. The fda [START_REF] Ramsay | fda: Functional Data Analysis[END_REF]) R (R Core Team 2020) package implements these methodologies and tools for dealing with such functional data.

In this paper we consider the case where the underlying stochastic model generating the data is a continuous-time stochastic process X = {X t , t ∈ T } such that for all t ∈ T , X t is a categorical random variable rather than a real-valued one.

Let (Ω, A, P) be a probability space, S = {s 1 , . . . , s K }, K ≥ 2, be a set of K states and X = {X t ; X t : Ω -→ S, t ∈ T }

(1) be a family of categorical random variables indexed by T . Thus, for some ω ∈ Ω, a path of X, X(ω), is a sequence of states s i j = s i j (ω) and time points t i = t i (ω) of transitions from one state to another one: {(t 0 , s i 0 ), (t 1 , s i 1 ), (t 2 , s i 2 ), . . .}

(2)

where 0 = t 0 < t 1 < t 2 < . . . are the jump times in T and s i j ∈ S with i j ∈ {1, . . . , K}, ∀j ≥ 0. This path is read as follows. At time t 0 = 0, ω is in some state s i 0 ; at time t 1 , t 1 > t 0 , ω moves randomly from s i 0 to the state s i 1 ; then at time t 2 > t 1 it moves from the state s i 1 to state s i 2 and so on. If T is the interval of time [0, T ] for some T > 0, then the observation process stops when the time T is reached or some absorbing state is observed.

We call the sample paths of X given by sequences of type (2) categorical functional data generated by the process X.

Figure 1 presents the graphical representation of one observation of categorical functional random variable. The representation (a) is appropriate when no natural order relationship exists on the set of states S, whereas the representation in (b) supposes that there exists some order relationship (≺) on S: s 1 ≺ s 2 ≺ . . . ≺ s K .

Figure 1: Categorical functional data: graphical representation.

To the best of our knowledge, and quite surprisingly, there are no recent researches devoted to this type of functional data despite its ability to model real situations in different fields of applications: health and medicine (status of a patient over time), economy (status of the market), sociology (evolution of social status), and so on. As a start point of research on this topic we mention the works of [START_REF] Boumaza | Contribution à l'Étude Descriptive d'une Fonction Aléatoire Qualitative[END_REF], [START_REF] Deville | Analyse de Données Chronologiques Qualitatives : Comment Analyser des Calendriers ?[END_REF], [START_REF] Deville | Correspondence Analysis with an Extension towards Nominal Time Series[END_REF] and [START_REF] Saporta | Méthodes Exploratoires d'Analyse de Données Temporelles[END_REF]. These works are devoted to the extension of factorial techniques (canonical and multiple correspondence analysis) towards functional data. Applications of these techniques are presented in [START_REF] Heijden | Multiple Correspondence Analysis as a Tool for Quantification or Classification of Career Data[END_REF] for analysing career data and in [START_REF] Preda | Analyse Harmonique Qualitative des Processus Markoviens des Sauts Stationnaires[END_REF] for studying spectral properties of the transition probability matrix of a stationary Markovian jump process with continuous time. In [START_REF] Cardot | Estimating Finite Mixtures of Semi-Markov Chains: an Application to the Segmentation of Temporal Sensory Data[END_REF] the authors are interested to cluster paths of semi-Markov processes using mixtures with application to sensory data.

In this article we present categorical functional data analysis as an extension of the multiple correspondence analysis towards functional data and its implementation in the cfda R package.

The theoretical foundations of this work are given in [START_REF] Deville | Analyse de Données Chronologiques Qualitatives : Comment Analyser des Calendriers ?[END_REF] and are based on the concept of optimal encoding of the states of the process X with respect to maximum variance criterion among all encodings. In Section 2 we present the theoretical background of the optimal encoding methodology defining the principal components of the process X throughout the optimal encodings. The approximation of the optimal encodings of the states into a basis of functions and optimal representation of categorical functional data in lower dimensional spaces are detailed. The implementation of the optimal encodings is presented throughout the cfda R package in Section 3 where an application on real data set (care trajectories for patients diagnosed with severe infection) is performed in view of visualisation, descriptive statistics and clustering.

Categorical Functional Data Analysis

Introduced in [START_REF] Saporta | Méthodes Exploratoires d'Analyse de Données Temporelles[END_REF] and [START_REF] Deville | Analyse de Données Chronologiques Qualitatives : Comment Analyser des Calendriers ?[END_REF] under the name "analyse harmonique qualitative", multiple correspondence analysis is extended to categorical functional data. There are several ways to do this, we have chosen in this work to introduce it as a problem of finding the latent variables (principal components) that are the most related to the process X = {X t , t ∈ T }. Therefore, the principal components will allow to define optimal encoding of states S = {s 1 , . . . , s K }.

Without loss of generality, let suppose that T = [0, T ], with T > 0. For x, y ∈ S and ∀t ∈ [0, T ], let denote by:

•

1 x t =      1 if X t = x, 0 otherwise, • p x (t) = P(X t = x) and p x,y (t, s) = P(X t = x, X s = y).
The general hypotheses considered in that framework are:

H 1 : the process X is continuous in probability, lim h→0 P(X t+h = X t ) = 0
and H 2 : for each time t ∈ [0, T ] (except possibly a finite discrete set of time points), any state has a strictly positive probability to occur:

p x (t) = 0, ∀x ∈ S, ∀t ∈ [0, T ].

The Principal Components

Let L 2 (Ω) be the space of real random variables with finite second moment and, for some t ∈ [0, T ], L(X t ) be the linear space spanned by X t . Then, the conditional expectation operator associated to X t ,

E t : L 2 (Ω) → L(X t ), z ∈ L 2 (Ω), z -→ E t (z) = x∈S E(z|X t = x)1 x t ,
is also the orthogonal projector on the space linearly spanned by the set of indicator random variables {1 x t , x ∈ S}. Notice that E t is self-adjoint, idempotent and of rank K.

For z ∈ L 2 (Ω) and t ∈ [0, T ], the coefficient η 2 (z; X t ) = VAR(E t (z)) VAR(z)
is a measure of the correlation between z and the variable X t . The empirical version of η 2 is known as Wilks' Lambda statistics, well known in multivariate ANOVA [START_REF] Nath | A New Statistic in the One Way Multivariate Analysis of Variance[END_REF].

Let us recall that if t 1 , t 2 , . . . , t p are p different time points in [0, T ], then the random variable z which maximizes

p i=1 η 2 (z; X t i ) (3)
defines the first principal component of the multiple correspondence analysis of the set of p categorical variables {X t 1 , X t 2 , . . . X tp } [START_REF] Escofier | Analyse Factorielle et Distances Répondant au Principe d'Équivalence Distributionnelle[END_REF]. By an iterative procedure, the principal components of higher order are defined as maximizing (3) under orthogonality conditions with respect to the principal components of lower order.

In [START_REF] Saporta | Méthodes Exploratoires d'Analyse de Données Temporelles[END_REF] and [START_REF] Deville | Analyse de Données Chronologiques Qualitatives : Comment Analyser des Calendriers ?[END_REF] the authors extend the multiple correspondence analysis to the process X = {X t , t ∈ [0, T ]} (seen as infinite set of categorical random variables). More specifically, the principal components are defined as the random variable z that maximizes the criterion

T 0 η 2 (z; X t )dt. ( 4
)
They show that under the hypotheses H 1 and H 2 , the variable z which maximizes (4) is the variable associated to the largest eigenvalue of the following (stochastic) eigenvalue problem:

T 0 E t (z)dt = λz. ( 5
)
The operator Q = T 0 E t dt is positive, hermitian, and compact. Therefore, Q has a countable set of positive eigenvalues and eigenvectors,

{(λ i , z i )} i≥1 such that λ 1 ≥ λ 2 ≥ . . . ≥ 0 and T 0 E t (z i )dt = λ i z i .
The variables {z i } i≥1 are called principal components of the process X. Notice that z = 1 (constant) is an eigenvector of Q associated to the largest eigenvalue λ max = T . It follows that the principal components {z i } i≥1 form a set of zero-mean uncorrelated random variables.

Moreover, we have that

i≥1 λ i = KT,
where K is the number of states. Thus, excluding the trivial eigenvalue λ max = T , the contribution of the i-th principal component z i to (4) is

Ctr(z i ) = λ i (K -1)T .

Optimal Encoding Functions

In order to solve (5), let denote by

ξ t = 1 λ E t (z), ∀t ∈ [0, T ]. ( 6
)
Under the hypotheses

H 1 and H 2 , for each t ∈ [0, T ] ξ t is X t -measurable i.e. E t ξ t = ξ t and {ξ t } t∈[0,T ] is a L 2 -continuous stochastic process. From (5) it follows that z = T 0 ξ t dt. ( 7
)
Taking the conditional expectation with respect to X t in (5), one obtains that the stochastic process {ξ t } t∈[0,T ] is eigenvector of the following (stochastic) eigenvalue problem posed in the space of L 2 -continuous stochastic processes:

T 0 K(t, s)ξ s ds = λξ t , ∀t ∈ [0, T ], (8)
with K(t, s) = E t E s , for all t, s ∈ [0, T ]. Recall that the spectral analysis of the kernel K(t, s) yields to the canonical analysis of X t and X s [START_REF] Hotelling | Relations Between Two Sets of Variates[END_REF].

It can be shown [START_REF] Saporta | Méthodes Exploratoires d'Analyse de Données Temporelles[END_REF]) that the eigenvalue problems ( 5) and ( 8) are equivalent in that sense that they have the same set of eigenvalues {λ i } i≥1 and there is an one-to-one correspondence between the principal components z i and the process ξ i = {ξ i,t , t ∈ [0, T ]}, ∀i ≥ 1. This correspondence is given by (6).

As in (5), the solution of ( 8) is unique up to a constant. To have unique eigenvectors, the usual constraint on {ξ t } t∈[0,T ] is that of total variance equals to one,

T 0 VAR(ξ t )dt = T 0 E(ξ 2 t )dt = 1. ( 9
)
The relation (9) implies

VAR(z) = E(z 2 ) = λ. ( 10
)
It follows from (6) that ξ t is X t -measurable for all t ∈ [0, T ] and one can write

ξ t = x∈S a x (t)1 x t , ( 11
)
where {a x } x∈S are deterministic functions on [0, T ] that we call optimal encoding functions. Introducing (11) in (8) one obtains the following eigenvalue equation,

T 0 y∈S p x,y (t, s)a y (s)ds = λa x (t)p x (t), ∀t ∈ [0, T ], ∀x ∈ S, (12)
where, p x (t) = P(X t = x) and p x,y (t, s) = P(X t = x, X s = y).

The integral system ( 12) is a more "classic" one than ( 5) and ( 8). Under the hypothesis H 1 and H 2 it admits the sequence of eigenvalues {λ i } i≥1 associated to the optimal encoding eigen-functions {a x i , x ∈ S} i≥1 . Notice that the constraint conditions in (9) are expressed now in terms of optimal encoding functions,

T 0 x∈S [a x (t)] 2 p x (t)dt = 1. ( 13
)
According to (7), for i ≥ 1, the i-th principal component z i is derived from the i-th optimal encoding functions {a x i } as

z i = T 0 x∈S a x i (t)1 x t dt, ∀i ≥ 1. (14)

Expansion Formulas and Dimension Reduction

As a summary of the previous section, the three equivalent eigen-problems stated in (5), ( 8) and ( 12) yield to the following elements of the analysis of X:

• the set of principal components {z i } i≥1 are zero-mean and uncorrelated:

-E(z i ) = 0, ∀i ≥ 1. -COV(z i , z j ) = λ i if i = j, 0 otherwise.
• the set of eigen-processes {ξ i = {ξ i,t , t ∈ [0, T ]}} i≥1 which generates the principal components by (7),

z i = T 0 ξ i,t dt, i ≥ 1,
are zero-mean and of unit total variance.

• the optimal encoding functions,

{a x i = {a x i (t), t ∈ [0, T ]}} x∈S,i≥1
. They generate the eigen-processes ξ i by (11),

ξ i,t = x∈S a x i (t)1 x t , ∀t ∈ [0, T ].
They satisfy the normalization condition (13).

Expansion Formulas. As an analogy to the Karhunen-Loève expansion for the scalar processes [START_REF] Deville | Méthodes Statistiques et Numériques de l'Analyse Harmonique[END_REF], the following expansion formulas hold [START_REF] Saporta | Méthodes Exploratoires d'Analyse de Données Temporelles[END_REF]):

• for the process

X = {X t , t ∈ [0, T ]} throughout the indicators 1 x = {1 x t , t ∈ [0, T ]}:
1 x t = i≥1 z i a x i (t) 1 p x (t) , ∀x ∈ S. (15)
• for the bivariate joint probability function, p x,y = {p x,y (t, s), t, s ∈ [0, T ]}: applying the Mercer theorem [START_REF] Mercer | Functions of Positive and Negative Type and their Connection with the Theory of Integral Equations[END_REF] to the integral equation ( 12), one has the following expansion formula:

p x,y (t, s) = p x (t)p y (s) i≥1 λ i a x i (t)a y i (s), ∀t, s ∈ [0, T ], ∀x, y ∈ S. (16)
In particular, for x = y and s = t we obtain

p x (t) =    i≥1 λ i [a x i (t)] 2    -1 , ∀t ∈ [0, T ], ∀x ∈ S. ( 17
)
Dimension Reduction. Using only the q first terms in the right-side part of ( 15) , q ≥ 1, one obtains the best approximation of order q of X (viewed as a vector process X = {1 x , x ∈ S}) under the L 2 norm, among all the linear expansions of type

1 x t ≈ q i=1 z i a x i (t) 1 p x (t)
, ∀x ∈ S.

Thus, the q first principal components, {z 1 , . . . , z q }, q ≥ 1, allow for -graphical representation of sample paths of X in R q (especially for q = 2, one obtains a 2-D representation of categorical functional data),

-fit of clustering and regression models with X as explanatory variables.

Approximation of Optimal Encoding Functions: a Basis Expansion Approach

The eigenvalue equation in ( 12) provides the optimal encoding functions. For a two-state process, [START_REF] Saporta | Méthodes Exploratoires d'Analyse de Données Temporelles[END_REF] considers the birth-and-death process on [0, 1],

X t = 0, if t < θ, 1, if t ≥ θ, ( 18
)
where θ is a random variable uniformly distributed on [0, 1]. The authors provide in this case explicit formulas for the eigenvalues {λ i } i≥1 , the optimal encoding functions {a x i } i≥1 , x ∈ S, and the principal components {z i } i≥1 . In [START_REF] Preda | Analyse Harmonique Qualitative des Processus Markoviens des Sauts Stationnaires[END_REF], the author considers the case of stationary Markovian continuous time processes with reversible distribution. In this case, the system in ( 12) reduces to a system of linear second-order differential equations with constant coefficients.

In general, the solution of ( 12) is obtained by approximation. In his seminal work, [START_REF] Deville | Analyse de Données Chronologiques Qualitatives : Comment Analyser des Calendriers ?[END_REF] proposes to approximate the encoding functions {a x i } i≥1 , x ∈ S, into a basis of functions of dimension m, m ≥ 1. As in the classical framework of functional data [START_REF] Ramsay | Functional Data Analysis[END_REF]), the choice of m is a trade-off between complexity computation and precision of the approximation. In our simulation study (section 4) we show the influence of the choice of m on the approximation of optimal encodings.

Let {φ 1 , . . . , φ m }, φ i : [0, T ] → R, i = 1, . . . , m, be a basis of functions (Fourier, B-splines, monomial, etc.) and for each x ∈ S consider the approximation:

a x (t) ≈ α (x,1) φ 1 (t) + α (x,2) φ 2 (t) + . . . + α (x,m) φ m (t), ∀t ∈ [0, T ], (19)
where α x = α (x,1) , α (x,2) , . . . , α (x,m) ∈ R m is the column vector of the expansion coefficients of a x into the basis {φ 1 , . . . , φ m }.

Plugging ( 19) in ( 12) and ( 13) one obtains the following classical eigen-problem:

Gα = λF α, ( 20
) under the constraint α F α = 1, (21)
where α ∈ R m×K is the column vector obtained by the concatenation of the vectors {α x } x∈S and G and F are square matrices of size mK × mK defined as follows:

• The matrix G is the covariance matrix of the random variables {V (x,i) , x ∈ S, i ∈ 1, . . . m}, defined as

V (x,i) = T 0 φ i (t)1 x t dt, ∀x ∈ S, (22)
G = G (x,i),(y,j) = COV V (x,i) , V (y,j) , x, y ∈ S, i, j = 1, . . . , m , ( 23
) G = . . . (y, 1) . . . (y, j) . . . (y, m) . . .                                   • • • • • • • • • • • • • • • • . . . . . . . . . . . . . . . . . . . . . . . . (x, 1) • • • • • • • • • • • • • • • . . . . . . . . . . . . . . . . . . . . . . . . (x, i) • • • • • • • COV V (x,i) , V (y,j) • • • • • • • . . . . . . . . . . . . . . . . . . . . . . . . (x, m) • • • • • • • • • • • • • • • . . . . .
F = F (x,i),(y,j) = E U (x,i),(y,j) , x, y ∈ S, i, j = 1, . . . , m , (24)
where U (x,i),(y,j) is the random variable

U (x,i),(y,j) = T 0 φ i (t)φ j (t)1 x t 1 y t dt =        T 0 φ i (t)φ j (t)1 x t dt if x = y, 0 otherwise. ( 25
)
Thus, F is a block diagonal matrix, each block being a square matrix of size m × m corresponding to each x in S, {U (x,i),(x,j) , i, j = 1, . . . , m}.

Example. Let observe that if 0 = t 0 < t 1 < . . . < t m = T is a sequence of time points in [0, T ] and for i = 1, . . . , m, one defines the basis (B-splines of order 1),

φ i (t) =      1 if t ∈ [t i-1 , t i [, 0 otherwise then the random variable V (x,i) represents the time spent in the state x in the interval [t i-1 , t i [. Since V (x,i) = U (x,i),(x,i) , then F is a diagonal matrix with elements E(V (x,i) ), x ∈ S, i = 1, . . . , m.

Estimation

Notice that the random variables V (x,i) and U (x,i),(y,j) , x, y ∈ S, i, j = 1, . . . , m, are computed from X = {X t , t ∈ [0, T ]} throughout the basis of functions {φ 1 , . . . , φ m }.

Thus, if {X 1 , . . . , X n } is a sample of n paths of X corresponding to a random sample (ω 1 , . . . , ω n ) ∈ Ω n , then the corresponding samples V (x,i) (ω) and U (x,i),(x,j) (ω), ω ∈ {ω 1 , . . . , ω n }, provide two classical data sets, V and U, as:

• the V data set with n rows and Km columns for the V 's random variables,

V = ω V (s 1 ,1) • • • V (s 1 ,m) • • • V (x,i) • • • V (s K ,m) ω 1 V (s 1 ,1) (ω 1 ) • • • V (s 1 ,m) (ω 1 ) • • • V (x,i) (ω 1 ) • • • V (s K ,m) (ω 1 ) . . . . . . . . . . . . . . . . . . . . . . . . ω n V (s 1 ,1) (ω n ) • • • V (s 1 ,m) (ω n ) • • • V (x,i) (ω n ) • • • V (s K ,m) (ω n )
• and the U dataset with n rows and Km 2 columns for the U 's random variables, respectively:

U = ω • • • U (x,i),(x,1) • • • U (x,i),(x,m) • • • U (s K ,m),(s K ,m) ω 1 • • • U (x,i),(x,1) (ω 1 ) • • • U (x,i),(x,m) (ω 1 ) • • • U (s K ,m),(s K ,m) (ω 1 ) . . . . . . . . . . . . . . . . . . . . . ω n • • • U (x,i),(x,1) (ω n ) • • • U (x,i),(x,m) (ω n ) • • • U (s K ,m),(s K ,m) (ω n )
Therefore, the matrices G and F are estimated from the sample {X 1 , . . . , X n } by the matrices G and F , the covariance matrix estimator of the random variables V 's and the mean estimator of the random variables U 's. For each i and j in {1, . . . , m} and x and y in S one has:

G (x,i),(y,j) = COV V (x,i) , V (y,j) = 1 n -1 n h=1 V (x,i) (ω h )V (y,j) (ω h ) -n V(x,i) V(y,j) and F (x,i),(y,j) =      Ū(x,i),(y,j) = 1 n n h=1 U (x,i),(y,j) (ω h ) if x = y, 0 otherwise.
An estimate of i-th eigen vector of (20), α (i) , i ≥ 1, is the i-th eigenvector α i of the eigenequation ( 26),

G α = λ F α, ( 26
) under the constraint α F α = 1. ( 27
)
Then, for each x ∈ S, the i-th encoding eigen-function a x i is estimated by

a x i = m j=1 α i,(x,j) φ j , i ≥ 1. ( 28
)
The estimates for the encoding functions allow to compute the principal components z i for each unit ω in the sample (ω 1 , . . . , ω n ), as

z i (ω) = T 0 x∈S a x i (t)1 x t (ω)dt = x∈S m j=1 α i,(x,j) V (x,j) (ω), i ≥ 1. ( 29
)
Notice that the variance of z i equals the i-th eigenvalue λ i of (26),

VAR( z i ) = λ i , i ≥ 1. ( 30
)
Confidence bounds. Bootstrapping from the V and U datasets, throughout (26) one obtains an estimate of the covariance matrix of α i denoted with Σ i . Therefore, for each t ∈ [0, T ], we have

VAR(a x i (t)) = φ(t) Σ (i,x) φ(t),
where φ(t) is the column vector φ(t) = (φ 1 (t), . . . , φ m (t)) and Σ (i,x) is the covariance matrix of α i,x = α i,(x,1) , . . . , α i,(x,m) . Notice that Σ (i,x) is a submatrix of Σ i .

Then, for a confidence level 1 -u, u ∈ [0, 1], a confidence interval for a x i (t) is obtained as

CI 1-u (a x i (t)) = a x i (t) ± ζ 1-u 2 VAR(a x i (t)),
where ζ 1-u 2 is the quantile of order 1 -u 2 of the standard normal distribution.

The cfda Package through Examples

The cfda R package provides functions to analyze categorical functional data allowing to compute basic statistics such as transition tables or visualisation, and compute the optimal encodings. It uses ggplot2 [START_REF] Wickham | ggplot2: Elegant Graphics for Data Analysis[END_REF]

Data

Real Dataset

The cfda package is illustrated with the care dataset. It contains 2929 care trajectories for patients diagnosed with a severe infection. Each month from the diagnosis of the infection, the follow-up of each patient is filled in using one of the following 4 states: "D", the patient has not a medical follow-up, "C", the patient has a medical follow-up but no treatment, "T", the patient has a medical follow-up with a treatment but the infection is not suppressed and "S", the patient has a medical follow-up with a treatment and the infection is suppressed.

The dataset can be loaded running:

R> data(care) R> head(care, 10)

id time state 3 0 D 3 5 D 9 0 D 9 1 D 13 0 D 13 7 D 15 0 D 15 4 T 15 7 C 15 8 D
The dataset is in a specific format required by every function in this package. Data must be provided as a data.frame with 3 columns named id, time and state. the id column contains the identifiers of the different individuals, the time column the different time of records of state changes and the state column containing the state that occurs at the corresponding time. For example, in the dataset above, four individuals (3, 9, 13 and 15) are visible. The individual 15 has an initial state (D) at time t = 0, it stays in this state until t = 4 months, at which he moves to a new state (T), and so on.

Note that within each individual, the time values must be ordered. Concerning id and state, the used format is quite versatile: character, factor or integer can be used.

Generate a Dataset

There are two functions for generating datasets: generate_Markov and generate_2State.

The generate_Markov function generates individuals following a Markov process where different time values are defined by an exponential law. The function is defined as below:

generate_Markov(n, K, P, lambda, pi0, Tmax, labels)

with n the number of individuals, K the number of state, P the transition matrix of size K*K between states containing the probabilities of changing from one state to another, lambda a vector of length K containing the parameters of the exponential law associated with each state, pi0 a vector of length K containing the probabilities of being in the different states at time t = 0, Tmax the maximal time of each individual, labels a vector containing the state names.

The generate_2State function generates individuals starting at time t = 0 with state 0 and with with a unique change into the state 1 at a time t ∼ U(0, 1). The only argument of generate_2State is n the number of individuals.

Visualize a Dataset

The summary_cfd function gives an overview of the dataset by printing information such as the number of individuals, the time range or the number of states, etc. All the printed information are returned in a list.

R> summary_cfd(care) One can note that all individuals have the same start time value but not the same end time value. This does not meet the constraints for performing a functional data analysis. A subsample of the dataset must be extracted before performing the encoding computation.

A sample of the individuals from the care dataset is plotted using the plotData function (see Figure 2). Each line corresponds to an individual of the dataset, the successive changes of states are represented by different colors. plotData(data, col, addId, addBorder, sort)

The plotData function takes in argument a formatted data.frame (data) and additional aesthetic parameters:

group a vector of length the number of individuals of data containing a variable describing the group of each individual. Individuals from different groups are displayed on different subplots. Individuals whose group is NA are ignored.

addId a boolean to add the individuals id on the y-axis.

addBorder a boolean to add the black border around each individual.

sort a boolean to sort individuals according to the duration in their first state.

col allows users to customize state colors by providing a vector of the same length as the number of state. col is a character (named) vector containing defined color names from R (e.g. c("red", "blue", "darkgreen")) or RGB colors (e.g. c("#E41A1C", "#377EB8", "#4DAF4A")).

nCol only when group is used, the number of columns used to display the different groups.

R> plotData(care[care$id <= 100, ]) R> plotData(care, addBorder = FALSE, addId = FALSE, sort = TRUE)

Extract a Dataset Meeting the Constraints

To compute the encodings, each individual must have the same start and end time. This is not the case in the care dataset. So, we select patients with a follow-up of at least 18 months and works from t = 0 to t = 18 months.

First, we compute the duration of each individual using the compute_duration function.

It computes the last time value minus the first time value for every individual. It returns a named vector, with the id as names, containing the duration. The results can be plotted using the hist function with the output of compute_duration as argument; it returns a ggplot object that can be modified.

R> duration <-compute_duration(care) R> head(duration) 3 9 13 15 18 21 5 1 7 32 18 5

R> hist(duration)

The resulting plot is displayed in Figure 3

Number of Jumps

The compute_number_jumps function counts the number of changes in each individual's state. It has two arguments: data, the dataset in the right format and countDuplicated, a binary value indicating if jumps in the same state must be ignored (FALSE) or not, the default is FALSE. An hist function is provided to plot the distribution of the number of jumps. For the care dataset, the number of jumps varies between 0 and 8 (cf. Figure 5), with most patients with less than 6 jumps.

R> hist(nJump)

The states between which these jumps occur are visible using the statetable function that counts the number of transitions between each pair of states. Transitions between identical states can be removed from the output We note that the state S (infection suppressed) is not an absorbing state indicating some patients have relapsed.

Probability to be in a State

The last interesting statistic is the probability to be in a state at a given time which is computed using the estimate_pt function. recorded state for any time greater than the last time entry; if TRUE, it considers that for any time greater than the last time entry, the records are missing; the default is FALSE. This is an important parameter when individuals have different ending time. This function returns a list of two elements: t, a vector containing the time values, pt, a matrix (with the states in rows and the time values in columns) containing the computed probabilities. A plot function is provided to display the results, the first parameter is the output of estimate_pt function, the second is ribbon. If ribbon = FALSE, the probability for a state is displayed with a line, if TRUE, with a ribbon (cf. Figure 6). In this figure, we note that the probability of not having a follow-up (D) decreases with the time. The probability to be cured (S) has an opposite trend.

R> proba <-estimate_pt(care2) R> proba $pt 0 1 2 3 4 5 ... D 0.991 0.653 0.596 0.566 0.555 0.552 ... C 0.008 0.202 0.180 0.166 0.156 0.134 ... T 0.000 0.099 0.171 0.203 0.159 0.128 ... S 0.001 0.046 0.053 0.065 0.131 0.185 ... $t [1] 0 1 2 3 4 5 6 ...

R> plot(proba, ribbon = TRUE)

Continuous-time Markov Chain

Assume the data come from a continuous-time Markov chain. Note S the set of states and X t ∈ S, the state at time t then for i, j ∈ S and t, s ≥ 0, then:

P (X s+t = j|X s = i, {X u = x u : 0 ≤ u < s}) = P (X s+t = j|X s = i) = P ij (t).
P ij (t) is the probability that the chain will be in state j in t time units, given it is in state i.

A continuous-time Markov chain is completely described by its transition matrix P = (P ij ) and λ i , i ∈ S the parameters of the exponentially distributed sojourn time in state i.

The estimate_Markov function estimates the transition matrix (P) and the λ parameter (lambda) associated with the mean sojourn time spent in each state through 1/lambda.

R> mark <-estimate_Markov(care2) R> mark The estimated process can be plotted as a diagram with the plot function displayed in Figure 7. Each node represents a state with its mean sojourn time. An arrow between two nodes defined a possible transition with its probability. R> plot(mark, main = "care: transition graph")

Optimal Encoding

The main contribution of cfda is the computation of an optimal encoding for categorical functional data performed by the compute_optimal_encoding function. The two main parameters are data, the dataset in the cfda format, and basisobj, a basisfd object created using the different create.*.basis functions from the fda package. It also performs bootstrap for computing a confidence interval on the computed encoding; associated parameters are computeCI, a logical indicating if bootstrap must be performed, nBootstrap, the number of bootstrap samples, and propBootstrap, the proportion of individuals used for each bootstrap sample. Other parameters are nCores the number of cores to use, verbose, if TRUE, some information are printed during the process. The compute_optimal_encoding function uses integrate (R Core Team 2020) to compute integrals, parameters for this function can be passed through ... in particular subdivisions, the number of subdivisions to estimate the integral.

R> set.seed(42)

R> basis <-create.bspline.basis(c(0, 18), nbasis = 10, norder = 4) R> fmca <-compute_optimal_encoding(care2, basis, nCores = 7) The main part of the computation time comes from the computation of V and U, these two steps are performed with parallel computation. Once these matrices computed, bootstrap is performed at a low computational cost. The output object of compute_optimal_encoding is a list containing:

######### Compute encoding #########
eigenvalues eigenvalues of the problem (26) alpha coefficients of the different encoding for each eigenvector (a list of matrices) (26) pc principal components for each eigenvector F F matrix (see equation ( 24))

V V matrix (see equation ( 22))

G covariance matrix of V (see equation ( 23)) basisobj basisobj parameter bootstrap encoding for each bootstrap sample varAlpha a list containing Σ (i,x) ∀i, ∀x ∈ S, covariance matrix of α i,x = α i,(x,1) , . . . , α i,(x,m)

Plot Functions

Three plot functions are associated with the compute_optimal_encoding function, the first argument of these functions is the output of compute_optimal_encoding. The first one, the plot function plots the encodings associated with a given eigenvector (harm parameter, by default, the encodings associated with the first eigenvector are plotted). If compute_optimal_encoding was run with parameter computeCI = TRUE, then the confidence interval can be added on the plot using the parameter addCI = TRUE. A subset of the states can be plotted by providing a vector with the state names to the states parameter.

The plotEigenvalues function plots the computed eigenvalues. It has two extra boolean parameters: cumulative, if TRUE, the cumulative sum of the eigenvalues is plotted and normalize, if TRUE, eigenvalues are normalized such that their sum is equal to 1. The last one is the plotComponent function that plots the individuals using the given components (comp parameter, a vector of length 2 containing the components' number). The other arguments are addNames that adds the individual's names on the plot and some parameters to adjust the position and size of these names (nudge_x, nudge_y and size).

The plots for the care dataset are shown in Figure 8 and 9 and are produced by the following code:

R> plotEigenvalues(fmca, cumulative = TRUE, normalize = TRUE) R> plotComponent(fmca, comp = c(1, 2), addNames = FALSE) R> plot(fmca) R> plot(fmca, addCI = TRUE)

Extract the Encoding Functions

The computed encoding functions can be extracted using the get_encoding function as an object of class fd (functional object from fda) using fdObject = TRUE or as a matrix using : chr [1:10] "bspl4.1" "bspl4.2" "bspl4.3" "bspl4.4" ... ..-attr(*, "class")= chr "basisfd" $ fdnames:List of 3 ..$ args: chr "time" ..$ reps: chr [1:4] "reps 1" "reps 2" "reps 3" "reps 4" ..$ funs: chr "values" -attr(*, "class")= chr "fd" We decide to keep 4 clusters regarding the heights of the tree. The resulting clusters can be displayed using the plotData function with the group argument.

R> class <-cutree(hc, k = 4) R> plotData(care2, group = class, addId = FALSE, addBorder = FALSE, + + sort = TRUE)

The different clusters are associated with the time spent in the different states after leaving the state "D" (cf. Figure 12). For example, the cluster number 1 corresponds to individuals that have spent most of their time (after "D") in the "C" state.

Simulation Study

Birth-and-death Process

Data are simulated under the simple model of birth-and-death process presented in [START_REF] Saporta | Méthodes Exploratoires d'Analyse de Données Temporelles[END_REF]. The process is defined on the interval time [0, 1] by

X t = 0, if t < θ, 1, if t ≥ θ, ( 31
)
where θ is a random variable uniformly distributed on [0, 1]. Table 1: True and estimated eigenvalues for the birth-and-death process. The estimated values are the mean over 100 samples. In brackets, the standard error is displayed. m=5 true n=50 n=100 n=200 n=500 1 0.5000 0.5117 (6.4e-4) 0.5013 (4.4e-4) 0.5000 (3.1e-4) 0.5009 (1.8e-4) 2 0.1667 0.1680 (3.3e-4) 0.1672 (2.0e-4) 0.1679 (1.5e-4) 0.1664 (0.9e-4) 3 0.0833 0.0824 (1.9e-4) 0.0835 (1.4e-4) 0.0834 (0.9e-4) 0.0835 (0.5e-4) 4 0.0500 0.0455 (1.3e-4) 0.0490 (0.9e-4) 0.0494 (0.6e-4) 0.0492 (0.4e-4) 5 0.0333 0.0184 (0.6e-4) 0.0205 (0.4e-4) 0.0211 (0.3e-4) 0.0215 (0.2e-4) m=10 true n=50 n=100 n=200 n=500 1 0.5000 0.5124 (6.4e-4) 0.5016 (4.3e-4) 0.5002 (3.1e-4) 0.5009 (1.8e-4) 2 0.1667 0.1692 (3.4e-4) 0.1677 (2.0e-4) 0.1682 (1.5e-4) 0.1665 (0.9e-4) 3 0.0833 0.0841 (2.0e-4) 0.0843 (1.4e-4) 0.0839 (0.9e-4) 0.0837 (0.5e-4) 4 0.0500 0.0486 (1.3e-4) 0.0510 (0.9e-4) 0.0508 (0.7e-4) 0.0501 (0.4e-4) 5 0.0333 0.0317 (0.8e-4) 0.0335 (0.6e-4) 0.0334 (0.4e-4) 0.0336 (0.2e-4) m=20 true n=50 n=100 n=200 n=500 1 0.5000 0.5128 (6.4e-4) 0.5018 (4.4e-4) 0.5003 (3.1e-4) 0.5010 (1.8e-4) 2 0.1667 0.1699 (3.4e-4) 0.1681 (2.0e-4) 0.1683 (1.5e-4) 0.1666 (0.9e-4) 3 0.0833 0.0849 (2.0e-4) 0.0847 (1.4e-4) 0.0841 (0.9e-4) 0.0837 (0.5e-4) 4 0.0500 0.0495 (1.3e-4) 0.0514 (0.9e-4) 0.0510 (0.6e-4) 0.0502 (0.4e-4) 5 0.0333 0.0327 (0.8e-4) 0.0340 (0.6e-4) 0.0337 (0.5e-4) 0.0337 (0.2e-4) • The approximation of optimal encoding functions in a basis of functions is based on the computation of random variables V (x,i) and U (x,i),(x,j) defined in ( 22) and ( 25) respectively. The computation of integrals involved in the definition of these random variables uses the inprod function of the fda R package which, at its turn, calls the function eval.fd. For n and K fixed, this step is the most computational in terms of time resources and it depends of the number of basis functions, m, considered for the approximation (19). As the computation is at done for every ω in (ω 1 , . . . , ω n ), parallel computation in performed.

• The F matrix ( 24) can be singular in some specific situations, namely when there exists an interval I ⊂ [0, T ] and some state x such that p x (t) = 0, ∀t ∈ I. In this case the hypothesis H 2 is not satisfied. For t ∈ I, the operator E t is degenerated, however, the eigenvalue equation ( 5) is still valid. From ( 12), the optimal encodings a x function is not defined for t ∈ I. From a computational point of view, if φ i is some element of the basis functions {φ 1 , . . . , φ m } with support in I then the random variables V (x, i) and U (x,i),(.,.) are zero-constant and therefore, the row and column corresponding to (x, i) in the F and G matrices are zero vectors. Thus, the element α (x,i) of the expansion coefficients vector α x is not defined. Dropping the rows and columns from F and G corresponding to (x, i) allows to solve the eigen-problem in (26). Notice that the constraints ( 27) are fulfilled.

The results in this paper were obtained using R 4.0.3 with the cfda 0.9.7 package. R itself and all packages used are available from the Comprehensive R Archive Network (CRAN) at https://CRAN.R-project.org/.

Summary and discussion

Categorical functional data is represented by paths of a continuous-time stochastic process with values in a finite set of states. Less popular than the real-valued functional data, that is yet another kind of infinite dimensional object. The analysis of categorical functional data is presented in this paper as an extension of the multiple correspondence analysis for the finite dimensional setting. principal components, optimal encoding functions and optimal representations are presented. A simulation study and a real data application illustrate the methodology implemented in the cfda R package.

Figure 2 :

 2 Figure2: A sample (left) and all individuals (right) from the care dataset plotted using the plotData function.

Figure 3 :

 3 Figure 3: Distribution of the duration of trajectories.

  Figure4displays the graphic for the care dataset. We note that people tends to stay longer without medical follow-up (D) than in the other states.

Figure 4 :

 4 Figure 4: Distribution of time spent per state.

Figure 5 :

 5 Figure 5: Distribution of number of jumps per individual.

Figure 6 :

 6 Figure 6: Probabilities to be in each state with regards to the time.

Figure 7 :

 7 Figure 7: Transition graph displayed using plot.Markov.

  *********************************************** DONE in 1.3s Run Time: 149.84s

Figure 8 :

 8 Figure 8: Plots generated by different graphical functions on the output of the compute_optimal_encoding function.

  Figure 9: Plots generated by the plot function on the output of the compute_optimal_encoding function.

Figure 10 :Figure

 10 Figure 10: Individuals with extreme negative value (min) and extreme positive value (max) on the component 1.

Figure

  Figure Content of the different clusters.

Figure 13 :Figure 14 :

 1314 Figure 13: True (solid black) and estimated encoding functions for state 0 of the birth-anddeath process. The estimated encoding functions are the mean of 100 samples.

  cumulative distribution function (b) Cumulative distribution function for z 2

Figure 15 :

 15 Figure 15: Empirical (resp. theoretical) cumulative distribution function for the first (z 1 ) and second (z 2 ) principal components of the birth-and-death process for n = 500 and m = 20.

  package to display graphics and the snow[START_REF] Tierney | snow: Simple Network of Workstations[END_REF] and foreach[START_REF] Microsoft | foreach: Provides Foreach Looping Construct[END_REF] packages for code parallelization.Other packages for analyzing sequences of categorical data exist, but not in a functional way. The TraMineR[START_REF] Gabadinho | TraMineR: Trajectory Miner: a Toolbox for Exploring and Rendering Sequences[END_REF]) R package provides functions to perform descriptive analysis, and distance functions between sequences to perform clustering analysis. The WeightedCluster[START_REF] Studer | WeightedCluster Library Manual: A Practical Guide to Creating Typologies of Trajectories in the Social Sciences with R[END_REF]) package relies on TraMineR and implements a clustering method and some cluster quality statistics. The msm[START_REF] Jackson | Multi-State Models for Panel Data: The msm Package for R[END_REF]) package estimates a continuous-time (hidden) Markov multi-state model. The R packages ClickClust[START_REF] Melnykov | ClickClust: An R Package for Model-Based Clustering of Categorical Sequences[END_REF]) and clickstream (Scholz 2016) model these sequences as discrete Markov chains to cluster them with mixture models or k-means.

  . Most trajectories last less than 40 months. To restrict individuals to a maximal time value of 18 months, we use the cut_data function that has two parameters: data and Tmax, the maximal time value. After applying this function, all individuals have Tmax as ending time.An interesting statistic is the time spent in each state per individual that can be computed with compute_time_spent function. It returns a matrix with n rows (number of individuals) and K columns (number of states) with the computed time. A plot function is provided to plot the distribution of each state.

	2 15	4	T
	3 15	7	C
	4 15	8	D
	5 15	15	C
	6 15	18	C
	7 18	0	D
	8 18	2	S
	9 18	18	S
	10 43	0	D
	3.2. Basic Statistics for Categorical Functional Data
	Time Spent in each State
	R> idToKeep <-names(duration[duration >= 18])
	R> care2 <-cut_data(care[care$id %in% idToKeep, ], 18)
	R> head(care2)
	id time state
	1 15	0	D

  table using removeDiagonal = TRUE.

	R> statetable(care2, removeDiagonal = TRUE)
	to
	from	D	C	T	S
	D	0 697 253 146
	C 271	0 346 97
	T 16 74	0 461
	S 16 91 31	0

Interpreting the Encoding Functions

First have a look at the plot of the encoding functions associated with the harmonic number 1 (cf. Figure 9a).

R> plot(fmca, harm = 1)

The lowest curve is for state "D", this indicates that individuals with a large negative value for component number 1 tend to spend more time in this state. Similarly, individuals with a large positive value tend to spend more time in the state "S". To check these statements, individuals with extreme values on the first component are plotted using the plotData function with the group parameter. A group variable is created with two different values: "min" for the individuals with the 5% lowest value, and "max" for the individuals with 5% highest value. R> minpc1 <-names(which(fmca$pc[,1] <= quantile(fmca$pc[,1], 0.05))) R> maxpc1 <-names(which(fmca$pc[,1] >= quantile(fmca$pc[,1], 0.95))) R> ids <-unique(care2$id) R> group <-factor(rep(NA, length(ids)), levels = c("min", "max")) R> group[ids %in% minpc1] = "min" R> group[ids %in% maxpc1] = "max" R> plotData(care2, group = group, addId = FALSE, addBorder = FALSE, + sort = TRUE) + + ggplot2::labs(title = "Extreme individuals on component 1")

The result is visible in Figure 10 and confirms our suppositions. Clearly, individuals in the "min" group spend 18 months in the "D" state (without medical follow-up) whereas individuals in the "max" group spend most of their time in the state "S" (infection suppressed).

Application to Clustering

The proposed method produces numerical encoding for categorical functional data. These encoding can be used for classical statistical methods such that regression or clustering. In the following, we perform a hierarchical clustering to find a structure in the care dataset.

The clustering is performed with the first principal components explaining at least 90% of the variance. The associated tree is displayed in Figure 11. R> nPc90 <-which(cumsum(prop.table(fmca$eigenvalues)) > 0.9)[1] R> hc <-hclust(dist(fmca$pc[, 1:nPc90]), method = "ward.D2") R> plot(hc, labels = FALSE) cfda: Categorical Functional Data Analysis In this case, K = 2, S = {s 1 = 0, s 2 = 1} and

and for t > s,

From ( 12) and ( 14), the authors provide explicit formulas for the eigenvalues {λ i } i≥1 , the optimal encoding functions {a x i } i≥1 , x ∈ S, and the principal components {z i } i≥1 as follows:

• the eigenvalues are given by:

• if P n (u) = 1 2 n n! d n du n (u 2 -1) n is the Legendre polynomials of order n, then the principal components z i corresponding to λ i , are given, up to a constant, by:

In particular, for i = 1, 2,

is uniformly distributed on -3 2 ; 3 2 and

Let observe that z 1 and z 2 are linearly uncorrelated but related by

showing some regularity of the 2-D representation of data throughout the plot {(z 1 (ω), z 2 (ω), ω ∈ Ω}.

• for i = 1, 2, the optimal encoding functions are given by:

and

Results

We simulate data from the above process with different number of trajectories (individuals), n = 50, 100, 200, 500 and a B-spline basis functions of order 4 with different number of basis functions m = 5, 10, 20 (equidistant knots). Simulation results are compared to the theoretical results presented in Section 4.1 to ensure the well-behaviour of the implemented method.

Eigenvalues

The first five eigenvalues for the different settings are compared to the eigenvalues from (32) in the table 1. The estimations are presented together with the associated standard errors in order to measure the impact of the choice of the sample size (n) and the dimension of the basis (m).

Encoding Functions

Figure 13 shows the mean over 100 samples of the first and second encoding functions for the state 0 for m = 5. The true encoding functions (34) and ( 35) are displayed in solid black line.

As for the eigenvalues, the best estimates are achieved with n = 200 and n = 500. The same conclusion holds for m = 10 and m = 20 but the number of basis does not seem to influence the accuracy.

Principal Components

In Figure 14, we check the relation between the first and second principal component (33). The theoretical equation is displayed in black whereas the computed principal components for a sample with n = 500 and m = 20 are in red. We note the closeness of the computed components with the theoretical equation.

In Figure 15, the cumulative distribution functions (cdf) for the two first principal components are displayed as well as their empiric equivalent for n = 500 and m = 20. As described above, z 1 follow an uniform distribution between -3 2 and 3 2 , the empiric cdf (in red) is closed to the theoretical one. The same representation is made for z 2 .

Computational details


Sarah Wassermann

Zied Pedro Casas

Ben Houidi

Alexis Huet

Michael Seufert

Nikolas Wehner

Joshua Schüler

Shengming Cai

Hao Shi

Jinchun Xu

Zied Ben Houidi

Tobias Hossfeld

Dario Rossi

Are you on Mobile or Desktop? On the Impact of End-User Device on Web QoE Inference from Encrypted Traffic

Keywords: Web QoE, Smartphone vs. Desktop, Network Monitoring, Machine Learning, SpeedIndex, Encrypted Traffic

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

The Web is one of the most relevant components of the Internet. The performance of the Web is highly relevant to the success of every online service, as it severely impacts the engagement and churn of users. The assessment of a web service as perceived by the end user can be realized through the corresponding Web Quality of Experience (QoE), which is very challenging to measure. Different from other services, such as video streaming or gaming, web browsing is a composite of numerous multimedia components and embedded services; loading a web page today requires tens of flows to get the various page resources located in diverse servers from different content providers. In this complex process, the network plays a key role impacting users' Web QoE, forcing Internet Service Providers (ISPs) to deploy effective means to monitor Web QoE as perceived by their customers.

The literature on web performance analysis presents a wide range of objective metrics capturing the performance of web pages, including metrics such as Page Load Time (PLT), SpeedIndex (SI), and Above the Fold Time (AFT). However, all these metrics require access to the application layer, which is hidden from the eyes of the ISP by the wide deployment of end-to-end network traffic encryption.

The analysis of Web QoE from purely in-network, encrypted traffic measurements is yet an under-explored problem; in fact, we have been the first recently addressing it, for the specific case of desktop web browsing [START_REF] Huet | Revealing QoE of Web Users from Encrypted Network Traffic[END_REF], [START_REF] Huet | Web Quality of Experience from Encrypted Packets[END_REF], using packet-level features as input. By using controlled page-load experiments, where network data is simultaneously collected with groundtruth Web-QoE metrics such as SI and AFT, we have shown the potential of using supervised Machine Learning (ML) to infer these metrics from features computed on the encrypted stream of packets. In this paper, we follow a similar approach to [START_REF] Huet | Revealing QoE of Web Users from Encrypted Network Traffic[END_REF], [START_REF] Huet | Web Quality of Experience from Encrypted Packets[END_REF], extending the analysis in multiple new directions:

1 -Multi-device models: we consider Web QoE not only for desktop devices, but include smartphone and tablet Web QoE, conceiving cross-device generalizable models. The lion's share of Internet-access devices today is smartphones, with nearly three quarters of the world population using exclusively their smartphones to access the Internet by 2025 [START_REF]Cisco Annual Internet Report[END_REF]. As we find in our results, a model trained only on desktop browsing data provides poor Web-QoE estimation performance when applied to smartphone and tablet measurements.

2 -Web-QoE Estimation: besides training regression models to estimate Web-QoE objective metrics -in particular SI -, we rely on real end-user data from previous studies [START_REF] Da Hora | Narrowing the gap between QoS metrics and Web QoE using Above-the-fold metrics[END_REF] to build Web-QoE classification models for subjective metrics -e.g., Mean Opinion (MOS) Scores.

3 -ML Benchmark: we present an extensive benchmark comparing the performance of different ML models for both estimation tasks -regression for SI inference, and classification for QoE estimation.

4 -Packet and Flow-level Models: we conceive models working either at the packet-level or at the flow-level; we show that the proposed flow-level models achieve highly similar performance to the packet-level ones, but using an order of magnitude less input features, thus showing strong potential for a practical monitoring solution.

The remainder of the paper is organized as follows. Section II overviews the related work on Web-QoE monitoring and analysis. Section III presents the overall modeling and data-generation approach, including a characterization of the produced datasets for this study. In Section IV, we introduce and evaluate the proposed packet-level ML models for two different tasks, including Web-QoE inference (SI) and classification of QoE ranges (excellent, good, poor); here we also show how models trained on single-device measurementse.g., desktop -, result in strong performance degradation when applied to other device types such as smartphones and tablets. In Section V we present a (per-task) unified, multi-device model, which achieves high inference/prediction performance with a strong generalization potential across desktop, smartphone, and tablet devices. In this section we also extend the multi-device models to operate at the flow-level, showing that the same performance can be achieved using much less, and easier to compute, input features. Finally, Section VI concludes this paper. As a conclusion, the proposed multi-device, flowlevel models lay the basis for a generalizable, multi-device, Web-QoE passive monitoring system.

II. RELATED WORK

Initial Web-QoE models were based on Page Load Times (PLT) [START_REF] Ibarrola | Web QoE evaluation in multi-agent networks: Validation of ITU-T G.1030[END_REF], [START_REF] Egger | Waiting Times in Quality of Experience for Web Based Services[END_REF], and are still broadly used in the practice to infer user satisfaction in web browsing, e.g, under ITU-T [START_REF]Estimating End-to-end Performance in IP Networks for Data Applications[END_REF]. However, research has demonstrated that PLT is a poor proxy to user perception of web-page loading times. Indeed, the actual web content visible to the user is usually displayed much earlier, as most web pages often stretch beyond the browser's viewport. Additional in-browser metrics have been accordingly devised to better suit the page display on the screen. An approach is the so-called Above the Fold Time (AFT), i.e., the time until the visible portion of a web page has been fully loaded, which has been also tested within traditional Web-QoE models [START_REF] Hoßfeld | Speed Index: Relating the Industrial Standard for User Perceived Web Performance to Web QoE[END_REF]. Newer Web-QoE metrics have been proposed recently, such as the SI, which considers the whole visual progress of the page loading, by processing a video capture of the screen. Besides single metric modeling, MLbased approaches have been explored [START_REF] Da Hora | Narrowing the gap between QoS metrics and Web QoE using Above-the-fold metrics[END_REF], [START_REF] Gao | Perceived Performance of Top Retail Webpages in the Wild: Insights from Large-scale Crowdsourcing of Above-the-fold QoE[END_REF] to model Web QoE from a combination of metrics.

Another direction in the literature proposes to understand how external components influence Web QoE. Prior work [START_REF] Sackl | The influence of network quality fluctuations on web qoe[END_REF], [START_REF] Sackl | Quantifying the impact of network bandwidth fluctuations and outages on web qoe[END_REF] has studied the impact of network-quality fluctuations and outages on user Web QoE. But besides network quality, other components influence Web QoE. They are linked to the specific web-page content -usability [START_REF] Varela | QoE in the Web: A dance of design and performance[END_REF], aesthetics [START_REF] Varela | Towards an understanding of visual appeal in website design[END_REF], etc. -as well as device type: desktop, smartphone, tablets [START_REF] Baraković | Survey of research on Quality of Experience modelling for web browsing[END_REF]. Important to our study, these papers show that smartphones and tablets have their own characteristics, not only regarding screen sizes, but also in terms of content rendering and web designs. Most of these papers focus on Web QoE in controlled, small-scale lab environments.

Others directly rely on in-browser metrics as a proxy to infer Web QoE, conducting large-scale active measurement campaigns. For example, the impact of multiple features such as transport protocols and network performance on PLT and AFT is studied in [START_REF] Saverimoutou | A 6-month Analysis of Factors Impacting Web Browsing Quality for QoE Prediction[END_REF], based on a set of 244 million measurements collected during 6 months for the top-10000 Alexa websites. Other papers also measured the impact of similar features on PLT and SI or AFT in different countries and different types of networks [START_REF] Asrese | Measuring Web Latency and Rendering Performance: Method, Tools, and Longitudinal Dataset[END_REF], including mobile ones [START_REF] Rajiullah | Web Experience in Mobile Networks: Lessons from Two Million Page Visits[END_REF]. Most of prior work has stayed at the application level, which is problematic for ISPs, who have no direct access to in-browser metrics. In recent years, TLS encryption has even narrowed the information that ISPs can collect from the network side, and previous approaches [START_REF] Ihm | Towards Understanding Modern Web Traffic[END_REF] based on deep packet inspection (DPI) and HTTP-traffic analysis are no longer applicable. Other papers [START_REF] Bocchi | Measuring the Quality of Experience of Web Users[END_REF], [START_REF] Trevisan | PAIN: A Passive Web performance indicator for ISPs[END_REF] developed correlated approximations to the SI metric, such as Byte/Object-Index [START_REF] Bocchi | Measuring the Quality of Experience of Web Users[END_REF] and Pain-Index [START_REF] Trevisan | PAIN: A Passive Web performance indicator for ISPs[END_REF], which can be computed from statistics of packet-and flow-level measurements, thus seamlessly operating with encrypted traffic. In recent work [START_REF] Huet | Revealing QoE of Web Users from Encrypted Network Traffic[END_REF], [START_REF] Huet | Web Quality of Experience from Encrypted Packets[END_REF], we took a step further to directly infer the SI metric, using ML techniques mapping network (encrypted) packetlevel traffic features to SI, in desktop devices.

When it comes to the specific case of Web-QoE monitoring in mobile devices, there have been multiple papers using machine learning [START_REF] Aggarwal | Prometheus: Toward Quality-of-Experience Estimation for Mobile Apps from Passive Network Measurements[END_REF]- [START_REF] Wassermann | Machine Learning Models for YouTube QoE and User Engagement Prediction in Smartphones[END_REF] or simple modeling approaches [START_REF] Nikravesh | QoE Inference and Improvement Without End-Host Control[END_REF] to map application-layer metrics [START_REF] Casas | Next to You: Monitoring Quality of Experience in Cellular Networks From the End-Devices[END_REF], [START_REF] Wassermann | Machine Learning Models for YouTube QoE and User Engagement Prediction in Smartphones[END_REF] or network-QoS metrics [START_REF] Aggarwal | Prometheus: Toward Quality-of-Experience Estimation for Mobile Apps from Passive Network Measurements[END_REF], [START_REF] Balachandran | Modeling Web Quality of Experience on Cellular Networks[END_REF], [START_REF] Nikravesh | QoE Inference and Improvement Without End-Host Control[END_REF] into QoE-related metrics. From these, two papers [START_REF] Aggarwal | Prometheus: Toward Quality-of-Experience Estimation for Mobile Apps from Passive Network Measurements[END_REF], [START_REF] Balachandran | Modeling Web Quality of Experience on Cellular Networks[END_REF] are the closest to our work, but both propose analysis approaches which are no longer applicable due to HTTP-traffic encryption [START_REF] Balachandran | Modeling Web Quality of Experience on Cellular Networks[END_REF], or do not address the specific problem of web browsing [START_REF] Aggarwal | Prometheus: Toward Quality-of-Experience Estimation for Mobile Apps from Passive Network Measurements[END_REF].

In this paper, we deal with an unexplored and so far neglected issue: understanding the impact of the end-user device type on the analysis of Web QoE from in-network traffic measurements. To the best of our knowledge, this is the first paper addressing this challenge.

III. WEB-QOE DATASETS & MODELING APPROACH

The proposed solution to the Web-QoE monitoring problem consists of training supervised ML models to map networktraffic features, extracted from the encrypted network-webpage traffic, into relevant Web-QoE metrics. The approach is data-driven, and thus needs datasets containing both the collected traffic traces -the input -and the targeted Web-QoE metric -the ground truth. The diagram presented in Figure 1 presents the workflow and the different stages of the solution.

To fully control the generation of such datasets, we built a measurement testbed based on multiple private instances of WebPageTest (WPT) -https://www.webpagetest.org/, the default, open-source web-performance-analysis tool used both in industry and academia. Figure 1, rows (1) and (2) describe this testbed and its usage. Different from previous studies [START_REF] Huet | Revealing QoE of Web Users from Encrypted Network Traffic[END_REF], [START_REF] Huet | Web Quality of Experience from Encrypted Packets[END_REF],

[15]- [START_REF] Rajiullah | Web Experience in Mobile Networks: Lessons from Two Million Page Visits[END_REF], [START_REF] Bocchi | Measuring the Quality of Experience of Web Users[END_REF], which have studied Web QoE exclusively for desktop browsers and desktop devices (or in some exceptional cases, emulating mobile devices), our measurement testbed consists of three different, non-emulated types of devices, including a smartphone device (Google Pixel 2 XL), a tablet (Google Pixel Slate), and a desktop computer (laptop), using WPT agents for Android and Linux. Chrome (last stable version) is used as browser. Instead of leveraging in-device WPT traffic-shaping capabilities, devices are connected to the open Internet through independent network emulators (emu), which allows for more realistic network-access-performance configurations in terms of bandwidth, latency, packet-loss rate, etc. This allows for heterogeneity in the generated measurements. Configurations used in the study include accessdownlink bandwidth up to 10 Mbps, packet-loss rates up to 10%, and RTTs up to 100 ms. Using WPT measurements, the platform extracts about 90 different KPIs and Web-QoE metrics, indicated as L7 Web-QoE logs in Figure 1, row

, such as PLT, SI, AFT, Byte Index [START_REF] Bocchi | Measuring the Quality of Experience of Web Users[END_REF], and Time to Interactive (TTI), as well as content characteristics of the visited pages. Network traffic is captured at an intermediate passive monitoring device (netmon) and stored as .pcap traces, from where model input features are extracted, indicated as L3 network-traffic features in Figure 1, row [START_REF] Huet | Web Quality of Experience from Encrypted Packets[END_REF]. For this study, we generated a per-device-type balanced dataset of more than 50.000 web page loading sessions (i.e., the loading of a single page), targeting the top 500 websites according to the Alexa top-sites list (https://www.alexa.com/topsites). The same pages are visited multiple times for each device type, using the same access-network setups. We do not consider the effect of caching, i.e., tests correspond to a first-view loading session. As we have recently shown in [START_REF] Huet | Revealing QoE of Web Users from Encrypted Network Traffic[END_REF], while caching has an impact on the performance of inference models, this impact is limited, and models generalize well across different protocol and caching settings [START_REF] Huet | Revealing QoE of Web Users from Encrypted Network Traffic[END_REF]. We focus on the inference of one particular Web-QoE metric, the SI, which is today one of the most accepted metrics reflecting Web QoE. Nevertheless, the methodology applies to any other similar Web-QoE metric. As shown in [START_REF] Bocchi | Measuring the Quality of Experience of Web Users[END_REF], measuring the SI is cumbersome in terms of computational resources and might introduce bias in the data capturing/processing, mainly due to the video capturing and analysis. This is particularly critical on smartphones and tablets, which are generally resource-constrained; therefore, instead of focusing on the SI metric, we collect the so-called RUM SpeedIndex (RUMSI) metric [START_REF] Meenan | Real User Monitoring SpeedIndex (RUMSI)[END_REF], which is a passive approximation to the SI, computed from the analysis of webpage resource timings.

To conclude, note that we assume that the measurement system takes as input network traffic from single web sessions. In an operational deployment in the wild -see row (3) in Figure 1 -, the traffic mix of concurrent web sessions must first be disentangled (step A) to then extract features from the traffic belonging to each web session (step B) and apply the trained model(s) (step C). This paper exclusively addresses steps B and C. Nevertheless, in case of concurrent web sessions, a classification methodology from the literature [START_REF] Trevisan | PAIN: A Passive Web performance indicator for ISPs[END_REF] could be applied to disentangle them. While the Web-traffic identification/disentangling (step A) is out of the scope of this paper, we have conceived multiple techniques addressing this problem.

A. Data Characterization

The list of top 500 Alexa pages is assorted in terms of contents, and as we show next, the type of device being used has a visible impact on web-page characteristics and timing performance. Figure 2(a) depicts the distribution of three relevant Web-QoE metrics, including the Time to First Paint (TTFP), which accounts for the time at which the first object is painted on the browser, the RUMSI, and the PLT. Note how the values are significantly higher for both smartphone and tablet devices as compared to desktop devices, pointing to a more complex rendering process in mobile devices. This is most probably linked to the specific hardware limitations of smartphones and tablets, as well as the particular characteristics of the OS and browser combination -native Chrome in Android. In addition, the way pages are optimized (most times dynamically) and rendered in mobile devices impacts loading times. Worse loading performance in mobile devices is a commonly known issue in practice; see page-speed statistics at https://backlinko.com/page-speed-stats. Interestingly, when comparing Android devices, TTFP values are almost identical for smartphone and tablet, RUMSI is slightly higher for tablet, while PLT is significantly higher for tablet. As we see next, this is most probably explained by the fact that web pages have more content to load in tablets. It is also interesting to note how PLT significantly overestimates the perceived loading time of web pages, represented by the (RUM)SI metric.

Figure 3 characterizes the 500 web pages per device type, in terms of (a) page size, (b) number of resources, (c) number of root domains, and (d) RUMSI to PLT ratio. The latter reflects the complexity of the web page in terms of visible content (SI) and full content downloading (PLT). As expected, web pages browsed in desktop devices are bigger than those browsed on smartphones or tablets, which are optimized for smaller screen sizes. The average page size is 2.7MB in desktop, 2.4MB in tablet, and 2.1MB in smartphone. Figures 3(b) and 3(c) further illustrate the richness and complexity of the web pages in terms of number of embedded contents and their location at different root domains, with more than 30% to 35% of the web pages consisting of more than 100 resources, and about 40% of the web pages fetching resources from more than 10 different root domains. The screen size probably plays a key role in terms of page characteristics, as the number of resources is higher for desktop, followed by tablet, and finally by smartphone. The RUMSI/PLT ratio shows not only how large the overestimation introduced by PLT is in terms of perceived page load times, but also how different this is for the different web pages. Indeed, about 10% of the pages have a ratio below 0.3 (the visible content loads way faster than the full content) and less than 5% of the pages have a ratio above 0.9 (the visible content basically corresponds to the full web page content). Conclusion: while most of the analyzed pages are very similar for every device type in terms of size and retrieved contents from external servers, differences can be significant for a small share of the pages. In terms of performance, loading times in Android devices (smartphone and tablet) are significantly higher than in desktop, a common trend observed in practice of web page speed analysis.

B. Subjective QoE Analysis

While the SI is a good objective metric reasonably capturing the Web QoE of real users [START_REF] Hoßfeld | Speed Index: Relating the Industrial Standard for User Perceived Web Performance to Web QoE[END_REF], we resort to previous subjective Web-QoE studies to better understand the expected QoE for the generated dataset. In particular, prior work [START_REF] Da Hora | Narrowing the gap between QoS metrics and Web QoE using Above-the-fold metrics[END_REF] conducted a subjective study where about 240 participants rated their browsing experience -loading of individual pages using a desktop browser -, according to a 5-level Absolute Category Rating (ACR) MOS score (bad QoE being 1 to excellent being 5). We rely on their publicly available dataset to identify QoErelated timing thresholds which could translate the RUMSI in our dataset to broad QoE classes. In Figure 2(c) we depict the relationship between SI and MOS scores obtained in that study. While the SI metric was not directly measured, additional metrics such as the Byte Index where computed,

C. Targets and Input Features

We realize the Web-QoE monitoring solution through two different prediction tasks: (i) inference of the RUMSI metric, which corresponds to a regression task, and (ii) prediction of the Web-QoE class {e,g,p}, which corresponds to a 3-class classification task. We use the same input features in both tasks, derived from the stream of encrypted packets. To define input features, we follow the rationale behind the computation of the SI metric itself, which considers the whole progress of the page loading. We define the Cumulative Bytes Downloaded features CBD(i) ∆T , as the (normalized) cumulative number of bytes downloaded from the first collected byte at time t 0 up to time t = t 0 + i × ∆T , with i = 1, . . . , m. The CBD features track the download progress of the page bytes, using a time resolution ∆T . Figure 4 depicts examples of CBD features for different network configurations, both for desktop and smartphone devices, using m = 100 and ∆T = 100 ms. Pages loading faster have a CBD loading curve rising sharper and arriving to full loading earlier.

For this study, we take m = 100 samples, and three different resolutions for the computation of features, using ∆T = 50 ms, 100 ms, and 500 ms, for a total of 300 CBD features. Using different resolutions helps capture different phenomena in the downloading progress, which potentially impact the SI, as well as allowing to track different page-load durations, in this case up to 5, 10, and 50 seconds, respectively. We consider n = 11 additional input session features, related to the complete pageloading session; these include: full session duration (first to last packet), downlink/uplink session duration (first to last packet in downlink/uplink direction), total number of packets downlink/uplink/full, total number of bytes downlink/uplink/full, and session mean throughput downlink/uplink. While these are mostly packet-level features, we extend in Section V-A the study to the implementation of flow-level features, realizing highly similar results. Figure 5 depicts the linear correlation between these input features and the RUMSI metric, for different device types. Correlation values are rather high for all devices, with stronger correlations observed for CBD features between 5 seconds and 10 seconds, as well as for session-duration features. Figure 6 shows correlation values for the QoE classification problem, for all devices together; as expected, based on the considered time-thresholds, higher correlations are observed between 2 and 5 seconds.

IV. DESKTOP MODELS' LACK OF GENERALIZATION

We now focus on the Web-QoE estimation tasks. As a reference for performance evaluation, we begin by training and evaluating different ML models for the specific case of desktop measurements. Recall that desktop measurements represent the most common data source so far used in the Web-QoE literature [START_REF] Huet | Revealing QoE of Web Users from Encrypted Network Traffic[END_REF], [START_REF] Huet | Web Quality of Experience from Encrypted Packets[END_REF], [START_REF] Saverimoutou | A 6-month Analysis of Factors Impacting Web Browsing Quality for QoE Prediction[END_REF]- [START_REF] Rajiullah | Web Experience in Mobile Networks: Lessons from Two Million Page Visits[END_REF], [START_REF] Bocchi | Measuring the Quality of Experience of Web Users[END_REF]. We then apply the trained models to both smartphone and tablet data, to highlight the cross-device lack of generalization and poor performance realized by single-device models in such multidevice scenarios. As a general note on the evaluations in this paper, all performance results correspond to 5-fold cross validation. bagging trees, and boosting including XGBoost [START_REF] Chen | XGBoost: A Scalable Tree Boosting System[END_REF]. The list is completed by a plain Bayesian approach, and by the k nearest neighbors algorithm (kNN). We assess performance using 3 performance metrics for regression problems, including the absolute error (AE), the relative error (RE), and the linear correlation (PLCC). We take both mean (M) and median (m) values for the error metrics, to filter out significantly large errors. Figure 7(a) additionally depicts the distribution of the prediction errors. RF100 attains the best inference performance, with a median absolute error of 249 ms, and a median relative error of 15%. Absolute prediction errors are below 500 ms for more than 70% of the sessions. More than 85% of the RUMSI values are inferred with an error below one second. Similar performance is realized by smaller ensembles, e.g., RF10, ET10, and bagging, using 10 instead of 100 trees. Given the training-speed improvements obtained with the ET10 model, we take it as the underlying prediction model in subsequent evaluations.

A. RUMSI Inference on Desktop

B. QoE Classification on Desktop

Table II reports the classification benchmark results obtained for desktop. Again, RF100 provides the best results, with an overall accuracy (ACC) close to 87%. Recall (R) and precision (P) are above 90% for the excellent QoE class prediction, but good and poor performance classes tend to be confused by the predictor. Nevertheless, recall and precision are close to 80% for these QoE classes.

C. Lack of Generalization for Mobile Devices

Now that we have built the models for desktop, a natural question is: how good would these models perform in data collected from other device types? This is critical in practice, as a significant, and ever growing, share of web-browsing activity comes from mobile devices. Figure 7(b) depicts the distribution of inference errors per device type, using the ET10 model, trained exclusively on desktop data. Table III summarizes the corresponding performance metrics. There is a strong inference-performance degradation when applying the desktop model to both smartphone and tablet data. Median absolute errors almost triple as compared to desktop performance. The distribution of errors shows that the desktop model tends to underestimate the RUMSI metric when applied to other devices, not present at training time. One could argue that RUMSI on smartphone and tablet is significantly higher than on desktop (cf. Figure 2), thus the performance degradation could be linked exclusively to this mismatch. However, performance degradation is also significant when only considering smaller RUMSI values, with median errors more than doubling for the example case of smartphone -from 260 ms to 592 ms -, testing only for RUMSI below 5 seconds. Figure 8 shows how the aforementioned cross-device training and validation issues also hold when considering different device types, using the RUMSI inference as example. The figure reports the usual AE, RE and PLCC metrics arranged as a training/testing matrix, where rows correspond to the devicetype data used for training, and columns to the device-type data used for testing. While specialization improves inference performance -the matrix diagonal -, training a model on measurements from a particular device type and applying the resulting model on measurements from a different device type results in poor inference performance, for all devicetype combinations. Note also how the cross-device lack of generalization applies to mobile devices, which are closer in terms of characteristics (cf. Figure 2 and Figure 3); while the performance degradation is lower when considering cross-data from smartphone/tablet devices, it is still non-negligible.

Performance degradation is also significant for the QoE classification problem. Figure 9 reports the classification performance per device type, again using the ET10 model, trained on desktop data. Overall classification accuracy drops from 84% on desktop to 67% on smartphone and tablet. Recall for excellent QoE degrades only slightly, but strongly for the other classes, and most importantly, precision for excellent QoE also drops strongly, meaning that the model cannot correctly track the classification problem. Conclusion: RUMSI inference and QoE prediction can be properly realized using CBD and session-based features, extracted directly from the stream of encrypted bytes. However, models so far proposed in the literature for single device types [START_REF] Huet | Revealing QoE of Web Users from Encrypted Network Traffic[END_REF], [START_REF] Huet | Web Quality of Experience from Encrypted Packets[END_REF], [START_REF] Saverimoutou | A 6-month Analysis of Factors Impacting Web Browsing Quality for QoE Prediction[END_REF]- [START_REF] Rajiullah | Web Experience in Mobile Networks: Lessons from Two Million Page Visits[END_REF], [START_REF] Bocchi | Measuring the Quality of Experience of Web Users[END_REF] might not perform properly in the wild, where other devices than desktop machines are used for web browsing.

V. MULTI-DEVICE (FLOW) WEB QOE MODELS

Having shown the lack of generalization and the crossdevice issues introduced by per-device models, we take the most natural step to conceive multi-device Web-QoE models. Possible approaches include the usage of stacking/ensembles of specialized models [START_REF] Casas | GML Learning, A Generic Machine Learning Model for Network Measurements Analysis[END_REF], or the inclusion of a pre-processing device-type classification task, preceding the main inference/prediction task. However, the simplest approach, exposing models to all devices data at training time, already provides high accuracy and generalizes well across devices, which has high practical appeal as it simplifies deployment.

Considering multi-device RUMSI inference first, Table IV and Figure 10(a) summarize the performance attained by a single ET10 model, trained on all-device data. Compared to per-device specialized models (cf. the matrix diagonal in Figure 8), there is a marginal degradation for the corresponding multi-device model, and mainly observed for desktop, with an error increase close to 10%. Still, performance is consistent across the three device types, with an overall median absolute error of 327 ms, and a relative error of 14%. Overall, the generalization capabilities of the multi-device model outweigh the accuracy of the specialized models, making it a preferred choice for Web-QoE monitoring in operational deployments.

Considering multi-device QoE classification next, Figure 10(b) reports the performance obtained with a single ET10 model, trained on all-device data: again, a slight performance degradation compared to the specialized desktop model (cf. Table II), but yields significant gain in terms of generalization to mobile devices (cf. Figure 9). The overall model accuracy is 82.2%, with recall and precision values for excellent{e}, good{g}, and poor{p} QoE of about {89%, 73%, 80%} and {86%, 74%, 82%}, respectively.

Conclusion: multi-device models significantly improve generalization of the Web-QoE inference across different devices, with only slight under-performance as compared to specialized models. As such, multi-device models provide simple, more accurate, and more reliable monitoring capabilities in realistic web-browsing scenarios.

A. Multi-device, Flow-level Models

In the last part of the study we focus on improving the practical application of the proposed Web-QoE inference models. In particular, we explore the definition of new input features at the flow level, which could be easier to compute than the proposed packet-level features so far considered. We define a set of 21 flow-level features, using similar notions to the ones which guided the packet-level features. These include: (i) the total number of flows (all, downlink, uplink), (ii) the min/mean/median/max flow duration in downlink, (iii) the min/mean/median/max flow size in downlink, (iv) the min/mean/median/max flow byte-index in downlink, (v) the mean/median in-flow, average intra-packets time (MDT) in downlink, (vi) the mean/median/max flow throughput in downlink, and (vii) the Flow-Index (FI). The flow byte-index uses the standard definition of Byte Index (BI) [START_REF] Bocchi | Measuring the Quality of Experience of Web Users[END_REF], but considering only the packets belonging to a specific flow. The FI feature represents an extension to the BI, but using flow size and flow ending time instead of packet size and time. Figure 11 depicts the basic notions behind the calculation of the FI. Both the FI and BI are integrallike metrics, similar to the definition of the SpeedIndex. We refer the reader to [START_REF] Bocchi | Measuring the Quality of Experience of Web Users[END_REF] for a comprehensive definition of the BI and the concepts of integral metrics. Flow features are complemented by the 11 session-level features, previously defined in Section III-C, adding to a total of 32 input features, computable at the flow level. For the sake of completeness, Figure 12 reports the linear correlation between these flow and session input features and the RUMSI metric. Features related to flow duration, flow BI, FI, and session duration are the ones showing the highest correlation to the RUMSI. Figure 13 and Table V report and summarize the RUMSI inference performance achieved by a multi-device model, using as input the new set of 32 features. Results are comparable, and even slightly better for some device types, than those achieved by using packet-level features (cf. Table IV), with the paramount advantage of using an order of magnitude less features, and at an easier-to-compute and more scalable level.

VI. CONCLUDING REMARKS

We have tackled the problem of Web-QoE monitoring from the ISP perspective, relying on in-network, passive measurements. Empirical evaluations on a large, multi-device, heterogeneous corpus of Web-QoE measurements for the most popular websites shows that the proposed solution can infer the (RUM)SI as well as estimate Web-QoE ranges from innetwork traffic measurements with high accuracy. At the same time, the device type introduces a strong bias in the capabilities of Web-QoE inference models, causing models trained for single device types to badly generalize to other devices. We showed that this cross-device lack of generalization can be solved by properly training on data coming from a multitude of devices. Our findings raise awareness of the fact that models for Web-QoE monitoring must be exposed to multi-device measurements to achieve proper inference and prediction performance in real network deployments, something generally neglected in the literature. To the best of our knowledge, we are the first to unveil the strong impact that the device type has on such models. The definition of novel flow-level features which also realize highly accurate predictions is also a keycontribution of the study.

Fig. 1 :

 1 Fig. 1: Diagram and workflow of the proposed solution.

Fig. 2 :

 2 Fig. 2: Distribution of (a) TTFP, RUMSI, and PLT values, (b) QoE classes per device type, based on (c) real-user MOS scores.

Fig. 3 :

 3 Fig.3: Data characterization, per device type, including page size, number of resources, number of root domains, and SI/PLT ratio. The latter reflects the complexity of the web page in terms of visible content (SI) and full content downloading (PLT).

Fig. 4 :

 4 Fig. 4: Examples of CBD features or loading curves, using ∆T = 100 ms, and different Access-RTT (RTTA) setups.

7 PLCCFig. 5 : 8 PLCCFig. 6 :

 7586 Fig. 5: Correlation between input features and RUMSI.

Fig. 7 :

 7 Fig. 7: RUMSI inference performance. Models are trained using exclusively desktop measurements.

Fig. 8 :Fig. 9 :

 89 Fig.8: Cross-device inference performance, using per-device ET10 as underlying model.

Fig. 10 :

 10 Fig. 10: (a) RUMSI inference and (b) QoE classification performance, ET10 model training on all-device data.

▪FlowFIFig. 11 :Fig. 12 :

 1112 Fig. 11: Flow-level features -flow index. ▪ Linear correlation to RUMSI ▪ Most correlated features ▪ mean/median/max flow duration ▪ mean/median/max flow BI ▪ mean MDT & FI ▪ session duration Flow-level & Session-level Features' Correlation

Fig. 13 :

 13 Fig. 13: Multi-device RUMSI inference performance -flowlevel features.

TABLE I

 I

	: Benchmarking of different ML models for RUMSI
	inference, for desktop.
	which can be used as a good proxy to the real SI [19]. Based
	on these subjective QoE results, we define 3 Web QoE classes:
	(e)xcellent -MOS ≥ 4 -, (g)ood -3 ≤ MOS < 4 -, and
	(p)oor -MOS < 3 -, resulting in SI thresholds of 2 and
	4 seconds. Interestingly, the SI thresholds recommended in
	the industry as target for excellent Web performance vary
	between 1 second (desktop) and 3 seconds (mobile), which
	are in line with the proposed higher QoE class threshold of 2
	seconds. It is important to note that our thresholds are derived
	for the case of browsing on desktop devices, and one would
	expect higher thresholds for Web QoE in mobile devices.
	Nevertheless, for this study, we assume the same thresholds
	apply to the three device types. Figure 2(b) shows that about
	60%/40% of the loading sessions correspond to excellent
	QoE in desktop/mobile (smartphone and tablet) respectively,
	25%/30% to good QoE, and the remaining 15%/30% result in
	poor QoE.

  Table I reports the RUMSI inference performance achieved by nine different ML models, most of them based on decision trees. The tested models include single decision tree (DT), multiple types of ensembles using different numbers of trees, such as extremely randomized trees (ET), random forest (RF),

	model				desktop
		ACC	R{e}	R{g}	R{p}	P{e}	P{g}	P{p}
	DT	80.3	88.8	66.1	72.8	88.4	66.2	73.8
	ET10	84.4	93.1	70.5	75.6	89.6	73.7	81.7
	RF10	84.6	93.1	71.6	74.6	90.1	73.1	82.1
	RF100	86.9	93.3	77.4	79.1	92.3	76.2	84.7
	Bagging	85.7	93.2	74.3	76.8	90.8	74.6	84.8
	Boosting	82.9	91.1	70.3	73.7	90.1	69.1	79.3
	Bayes	60.1	93.0	11.6	19.5	63.2	38.3	46.7
	kNN	74.9	87.3	55.1	62.1	81.8	59.1	72.0
	XGB	82.2	90.9	69.0	72.1	89.8	67.9	77.6

TABLE II :

 II Benchmarking of different ML models for Web-QoE prediction, on desktop. The three levels of QoE correspond to excellent{e}, good{g}, and poor{p} QoE.

	device	MAE-mAE (ms)	MRE-mRE (%)	PLCC
	desktop	598 -260	31 -16	0.879
	smartphone	1245 -721	41 -28	0.728
	tablet	1434 -724	44 -28	0.618
	smartphone	867 -592	43 -29	0.455
	RUMSI < 5s
	smartphone	2812 -2000	32 -27	0.667
	RUMSI > 5s

TABLE III :

 III Inference performance per device type. The ET10 model is trained using desktop data.

TABLE IV :

 IV Multi-device RUMSI inference performance.

TABLE V :

 V Multi-device, flow-level RUMSI inference.