Amay Cheam

Marc Fredette

Matthieu Marbac

Fabien Navarro

Translation-invariant functional clustering on COVID-19 deaths adjusted on population risk factors

Keywords: Functional data, Mixture models, Semiparametric models, Single-index regression, Wavelets

The COVID-19 pandemic has taken the world by storm with its high infection rate. Investigating its geographical disparities has paramount interest in order to gauge its relationships with political decisions, economic indicators, or mental health. This paper focuses on clustering the daily death rates reported in several regions of Europe and the United States over seventeen months. Several methods have been developed to cluster such functional data. However, these methods are not translation-invariant and thus cannot handle different times of arrivals of the disease, nor can they consider external covariates and so are unable to adjust for the population risk factors of each region. We propose a novel three steps clustering method to circumvent these issues. As a first step, feature extraction is performed by translation-invariant wavelet decomposition which permits to deal with the different onsets. As a second step, single-index regression is used to neutralize disparities caused by population risk factors. As a third step, a nonparametric mixture is fitted on the regression residuals to achieve the region clustering.

Introduction

In March of 2020, the World Health Organization (WHO) declared pandemic status for the novel coronavirus SARS-Cov-2, denoted COVID-19, indicating that it has reached a critical level of spread and severity worldwide. The global nature of COVID-19 pandemic has resulted in plenty of heterogeneity of the data, aggravated further by lack of prior knowledge or coordinated mitigation strategies which impeded research efforts. For instance, the assumption that the first occurrence emerged concurrently everywhere is improper. Additionally, the number of confirmed cases depends on the amount of tests that are being performed in a region. Hence, a region that has tested very few people can only report very few confirmed cases. Alternatively, the number of COVID-19 deaths is more systematically recorded: countries are asked to follow the 'cause of death' classifications from the WHO's International Classification of Diseases guidelines (World Health Organization, 2016). Though each country is responsible to provide their own guidance on how and when COVID-19 deaths should be recorded, this metric remains more reliable. Undoubtedly, the rapid propagation of this acute infectious respiratory disease has posed governmental challenges. Government responses to contain the virus's spread were multiple (social distancing, travel restrictions, lockdowns, etc.) and their efficiency needs to be investigated. To better understand this virus, it is profoundly useful to cluster regions similarly affected by This paper focuses on the investigation of spatial disparities of COVID-19 by considering the daily COVID-19 deaths recorded over seventeen months. Previous works consider this aim for analyzing a specific region [START_REF] Cioban | Spatial clustering behaviour of covid-19 conditioned by the development level: Case study for the administrative units in romania[END_REF] and propose to perform clustering by using geometrical methods (see [START_REF] Bullock | Mapping the landscape of artificial intelligence applications against covid-19[END_REF] and the reference inside). In this context, [START_REF] Bucci | Clustering spatio-temporal series of confirmed covid-19 deaths in europe[END_REF] propose to use a Bayesian nonparametric approach with a Gaussian mixture model to cluster regions in Europe, while [START_REF] Tang | Functional data analysis: An application to covid-19 data in the united states[END_REF] propose to use classical functional cluster analysis to cluster the continental states of the United States of America (USA). Alternatively, [START_REF] Chen | Clustering us states by time series of covid-19 new case counts with non-negative matrix factorization[END_REF] propose to use a non-negative matrix factorization (NMF) followed by a K-means clustering procedure applied on the coefficients of the NMF basis to cluster the continental states of the USA based on the new daily confirmed case counts. These approaches perform well when the study area is focused on a specific region. However, in this paper, we focus on clustering regions of the European Union (EU) and the USA. When considering regions of Europe and North America, a difficulty arises: COVID-19 outbreaks started at different times. The misalignment of the first occurrence between regions should not be neglected, whether between continents or within a country. To circumvent this issue, geographical disparities of the COVID-19 impact can be investigated by geometric clustering approaches based on dissimilarity matrices that are not sensitive to time shifting. In this context, the dynamic time warping (DTW) approach is the most popular method to obtain such dissimilarity matrix. Then, hierarchical clustering [START_REF] Park | Clustering county-wise covid-19 dynamics in north carolina[END_REF] or spectral clustering can be performed by considering this matrix [START_REF] Allem | A spectral clustering approach for the evolution of the covid-19 pandemic in the state of rio grande do sul, brazil[END_REF]. Note that [START_REF] Michael | Learning trends of covid-19 using semisupervised clustering[END_REF] the time shifting between the curves in a semi-supervised approach to learn different curve profile in UE and use these patterns to predict the evolution of the disease in the USA. Another problem to acknowledge is that the mortality occurs at different rates under different population risk factors [START_REF] Williamson | Factors associated with covid-19related death using opensafely[END_REF]. Hence the necessity to adjust these region-specific risk factors is intrinsic to allow regions to be compared fairly. Indeed, if these factors are neglected during clustering, then one observed a correlation between the estimated clusters and the risk factors [START_REF] Ramírez-Aldana | Spatial epidemiological study of the distribution, clustering, and risk factors associated with early covid-19 mortality in mexico[END_REF]. Furthermore, by adjusting the population risk factors, we are able to detect regions more susceptible to COVID-19 and perhaps identify the disparity factors between clusters. For instance, it allows for a retrospective assessment of the effectiveness and the quality of government responses, a concurrent analysis of the economic indicators, and a prospective perception of mental health during this unprecedented period.

Traditionally, clustering may be achieved through finite mixture models [START_REF] Mclachlan | Finite mixture models[END_REF]. When the family of distributions for each cluster is unknown, nonparametric mixtures can be considered to avoid unjustified parametric assumptions [START_REF] Chauveau | Semi-parametric estimation for conditional independence multivariate finite mixture models[END_REF]. A classical approach among these methods is to define the density of each mixture component as a product of univariate densities [START_REF] Hall | Nonparametric estimation of component distributions in a multivariate mixture[END_REF][START_REF] Kasahara | Non-parametric identification and estimation of the number of components in multivariate mixtures[END_REF]. The data we analyze are functional and thus raise the problem of data dimension [START_REF] Ferraty | Nonparametric functional data analysis[END_REF][START_REF] Ramsay | Applied functional data analysis: methods and case studies[END_REF]. To circumvent the curse of dimensionality, many model-based clustering approaches approximate the observed functions in some functional basis then perform clustering on the coefficients related to the basis (see the review of [START_REF] Jacques | Functional data clustering: a survey[END_REF] or [START_REF] Cheam | On the importance of similarity characteristics of curve clustering and its applications[END_REF]). For instance, [START_REF] Bouveyron | The discriminative functional mixture model for a comparative analysis of bike sharing systems[END_REF] proposed approximating the curves into a Fourier basis expansion coefficients, then perform clustering on the obtained coefficients with a Gaussian mixture. Alternatively, feature extraction could be accomplished via an orthogonal wavelet basis [START_REF] Antoniadis | Clustering functional data using wavelets[END_REF]. A further issue raised by this type of data is that of curve alignment. This has been addressed by previous works that do not tackle clustering [START_REF] Kneip | Statistical tools to analyze data representing a sample of curves[END_REF][START_REF] Wang | Alignment of curves by dynamic time warping[END_REF][START_REF] Ramsay | Curve registration[END_REF]. Recently, this issue has been considered for clustering in distance-based [START_REF] Paparrizos | k-shape: Efficient and accurate clustering of time series[END_REF] and model-based [START_REF] Chudova | Translation-invariant mixture models for curve clustering[END_REF][START_REF] Gaffney | Joint probabilistic curve clustering and alignment[END_REF][START_REF] Liu | Simultaneous curve registration and clustering for functional data[END_REF] frameworks.

In this article, we propose a novel three-step approach that circumvents the issues of clustering regions with respect to the COVID-19 dataset: the varying times of arrivals of the virus and the need to incorporate the population risk factors. This approach is named Clustering Regression residuals of Features given by Translation Invariant Wavelets (CRFTIW). The first step of CRFTIW consists of feature extraction using a multiscale approach based on translation-invariant (TI) wavelets [START_REF] Coifman | Translation-invariant de-noising[END_REF], which allows the shifted onsets of COVID-19 to be tackled by avoiding any pre-processing step for curve alignment (see [START_REF] Wang | Alignment of curves by dynamic time warping[END_REF] and the references cited in Jacques and Preda (2014, Section 2.3)). The objective to construct clusters that are invariant to time-shifts is somewhat different from conventional clustering in that it permits us to answer slightly different scientific questions about the data. Standard clustering (no time-shifts) will identify regions that peak at the same period, while TI clustering recovers regions that behave in similar patterns that unravel across time. The features are defined as the logarithm of the norm of the TI wavelet coefficients at each scale. The second step of CRFTIW integrates the population characteristics with a single-index regression of the features on the population risk factors. This approach has the benefits of the nonparametric regression but does not suffer from the curse of dimensionality. We show that the residuals of the regression preserve the cluster information. As the third step of CRFTIW, clustering of the regions is achieved by fitting a nonparametric mixture on the regression residuals. The only assumption made at this step is to define the density of each component as a product of univariates densities. The proposed approach has differences with the approach of [START_REF] Gaffney | Joint probabilistic curve clustering and alignment[END_REF] despite both methods considering curve translations. First, the scaling that we proposed depends on the covariates (i.e., the risk factors). Second, we use a wavelet approach that permits a greater reduction of the dimension. Finally, we consider a semi-parametric mixture that avoids the bias of the parametric mixtures observed when their parametric assumptions are violated.

The remainder of this article is organized as follows. Section 2 presents the data. Section 3 defines the new approach CRFTIW. Section 4 illustrates, with numerical experiments, the relevance of CRFTIW. In Section 5, we analyze the geographical disparities of the COVID-19 by investigating the regression coefficients, describing the clusters according to the diseases and illustrating the use of the clusters to investigate policy strategies. Section 6 presents some concluding remarks.

Description of the data

For this ongoing COVID-19 dataset, we consider n = 94 regions composed of the 50 states of the USA plus the District of Columbia and countries from Europe and America (Austria, Belgium, Bulgaria, Brazil, Canada, Switzerland, Chile, Colombia, Costa Rica, Cuba, Cyprus, Czechia, Germany, Denmark, Dominican Republic, Estonia, Spain, France, United Kingdom, Greece, Guyana, Croatia, Haiti, Hungary, Ireland, Iceland, Italy, Jamaica, Latvia, Netherlands, Norway, Panama, Peru, Poland, Portugal, Paraguay, Romania, Serbia, Sweden, Slovenia, Slovakia, El Salvador and Uruguay) whose data are available at Johns Hopkins' Github repository (file Policy.rds in [START_REF] Badr | Unified covid-19 dataset[END_REF]).

Our focus is on the curve W i = (W i(1) , . . . , W i(T ) ) recording the daily rate of the number of deaths per million people in each region i for a total of T = 512 days (between March 1st, 2020 to July 25, 2021, inclusively), where W i(t) denotes the death rate recorded for region i at time t. Data were extracted from the Center for Systems Science and Engineering at Johns Hopkins' Github repository (file Policy.rds in [START_REF] Badr | Unified covid-19 dataset[END_REF]) and a 7-day moving average has been performed due to the discrepancy of the data recorded by each region. For instance, this can account for days in the week where data may not be available, such as weekends. There are differences in time of arrival of the peak death rates between regions. This is illustrated by Figure 1 which shows that the first two COVID death peaks arrived first in Europe (see Austria and Italy), then in North America (see New Hampshire and Pennsylvania) and then in Latina America (Brazil and Costa Rica).

Early findings suggested that differences in COVID-19 disease prevalence and severity may be associated with certain risk factors [START_REF] Williamson | Factors associated with covid-19related death using opensafely[END_REF]. Thus, we consider two groups of risk factors. The first group contains three environmental risk factors (fine particulate matter PM2.5 concentration, nitrogen dioxide NO2 concentration and population density) and the second group contains eight medical risk factors (age-adjusted percent prevalence of adults with diagnosed diabetes, percent of obese adults, age-adjusted percent prevalence of adults who are current smokers, age-standardized percent prevalence of chronic obstructive pulmonary disease, age-standardized percent prevalence of cardiovascular disease, age-standardized percent prevalence of HIV/AIDS, percent of adults with hypertension and population proportion over 65 years old). All the risk factors were extracted from the Center for Systems Science and Engineering at Johns Hopkins' Github repository (file COVID-19 Static.rds in [START_REF] Badr | Unified covid-19 dataset[END_REF]) and have been scaled. For each group of risk factors, we perform a principal component analysis (PCA). For each region i, we store in X i the first principal components of both groups of risk factors. In the application described in Section 5, we show that it is relevant to consider the first principal component of the PCA based on the environmental risk factors and to consider the first two principal components of the PCA based on the environmental risk factors. Note that that we are able to consider the states of the USA because all the data (risk factors and daily rate of the number of deaths) are available at Johns Hopkins' Github repository, however, these data were not available for the different provinces of Canada or the states of Brazil.

Method

Outline of the three-step method

The daily COVID-19 death curves of the n regions W 1 , . . . , W n are supposed to independently arise from L different clusters. The cluster membership of region i is defined by the latent variable Z i = (Z i1 , . . . , Z iL ) where Z i = 1 if region i belongs to cluster and Z i = 0 otherwise. The model assumes that, conditionally on the cluster , each W i is defined as a product between a noisy version of δ i -lagged values of an unobserved curve u and the effect of the population risk characteristics µ(X i ) > 0, where X i denotes population risk factors of region i and where we set E[µ(X i )] = 1, for identifiability reasons. The deterministic functions u do not depend on the covariates X i . Moreover the noises ε i and the covariates X i are independent. Thus, given the cluster membership and the population risk factors of region i, we have

W i = L =1 z i µ(x i )(u (δ i ) + ε (δ i ) i ), (1)
where u (δ i ) and ε

(δ i ) i
are δ i -lagged versions of u and ε i , and the distribution of each ε i follows a centered distribution having a finite variance defined by the density

f (i.e., E f [ε i ] = 0 and E f [ε 2 i ] < ∞). Thus, noting that f (w i | X i = x i ) = L =1 f (w i | X i = x i , Z i = 1), f (w i | X i = x i , Z i = 1) = P(Z i = 1 | X i = x i )f (ε i | X i = x i , Z i = 1) = f ( w i µ(x i ) -u (δ i ) )
and, by assumption

P(Z i = 1 | X i = x i ) = P(Z i = 1) = π , the conditional distribution of W i given X i = x i is defined by the density f (w i | x i ) = L =1 π f w i µ(x i ) -u (δ i ) , (2)
where π > 0 is the proportion of cluster with L =1 π = 1. Note that the model defined by ( 2) is related to the translation invariant model proposed by [START_REF] Casa | Co-clustering of timedependent data via the shape invariant model[END_REF] for co-clustering of multivariate functional data. However, in [START_REF] Casa | Co-clustering of timedependent data via the shape invariant model[END_REF] the time of translation δ i is random variable that is relevant for clustering (i.e., the conditional distributions of δ i given the cluster memberships are not equal for all the clusters) while this variable is supposed to be irrelevant for clustering in (2). In addition, in [START_REF] Casa | Co-clustering of timedependent data via the shape invariant model[END_REF] the noises are supposed to be independent and identically distributed within a cluster and thus their distributions are not impacted by the time translation. In (2), the noises are not supposed to independent and identically distributed within components and thus their distributions can be impacted by the time translation.

Despite the model defined by (2) permitting a clustering of the regions based on the daily COVID-19 death curves with respect to the population risk factors, the estimation the multivariate density f is highly complex. Thus, we achieve the clustering with the following three-step approach:

1. Perform feature extraction of the daily COVID-19 death curves W i to obtain Y i ∈ R J+1 using TI wavelets (see Section 3.2).

2. Fit single-index regressions of the features Y i on the population risk factors X i and consider the residuals ξi ∈ R J+1 (see Section 3.3).

3. Use the nonparametric mixture to cluster the regions based on the residuals ξi (see Section 3.4).

This approach is relevant since the specific feature extraction reduces the dimension, permits us to deal with lagged values and keeps the main cluster information. Moreover, the singleindex regression keeps the cluster information of the features, allows for adjustment on population risk factors, and provides meaningful parameters used for detecting protective or compounding effects of the population characteristics and odd ratios.

Feature extraction and time misalignment

A wavelet basis is a set of functions obtained as translations and dilatations of two specific functions: a scaling function denoted by φ and a mother wavelet denoted by ψ. For the purpose of this paper, we use Daubechies wavelets and in particular the Symlet family.

We present the essentials below, more details can be found in [START_REF] Daubechies | Ten lectures on wavelets[END_REF] or [START_REF] Mallat | A wavelet tour of signal processing: the sparse way[END_REF]. Such wavelets are optimal in the sense that they have minimal support for a given number of null moments. For the simulations, we consider in particular the Symlet family (or Daubechies least asymmetric) which is a modified version, with more symmetry than the classical Daubechies wavelets. We present the essentials below, more details can be found in [START_REF] Daubechies | Ten lectures on wavelets[END_REF] or [START_REF] Mallat | A wavelet tour of signal processing: the sparse way[END_REF].

The scaling function φ and wavelet function ψ of an orthonormal Daubechies wavelet system satisfy the scaling equations

φ(t) = √ 2 k∈Z h k φ(2t -k), ψ(t) = √ 2 k∈Z g k φ(2t -k),
with scaling and wavelet filters h = {h k : k ∈ Z} and g = {g k : k ∈ Z}, with g k = (-1) k h 1-k . For any j ≥ 0, we set Λ j = {0, . . . , 2 j -1} and, for k ∈ Λ j ,

φ j,k (t) = 2 j/2 φ(2 j t -k), ψ j,k (t) = 2 j/2 ψ(2 j t -k).
Following the methodology of [START_REF] Cohen | Wavelets on the interval and fast wavelet transforms[END_REF], there exists an integer τ such that, for any integer j 0 ≥ τ , the collection of functions

S = {φ j 0 ,k , k ∈ Λ j 0 ; ψ j,k ; j ∈ N -{0, . . . , j 0 -1}, k ∈ Λ j }
forms an orthonormal basis of L 2 ([0, 1]). Therefore, for any integer j 0 ≥ τ and W ∈ L 2 ([0, 1]), we have the following wavelet expansion

W (t) = k∈Λ j 0 α j 0 ,k φ j 0 ,k (t) + ∞ j=j 0 k∈Λ j β j,k ψ j,k (t), t ∈ [0, 1], (3)
where α j,k and β j,k are the scaling and wavelet (or details) coefficients of f at scale j and position k defined as

α j 0 ,k = 1 0 f (t)φ j 0 ,k (t)dt, β j,k = 1 0 f (t)ψ j,k (t)dt.
The first term of this development can be seen as an approximation term of the function at level j 0 and the second as a detail term that characterizes the approximation error. The level of approximation j 0 has relatively no influence on the performance of the approximation (or estimation) and will subsequently be set to 0. The decomposition of the observations in a given wavelet basis is defined by

W i (t) = α i,0,0 φ 0,0 (t) + J-1 j=0 2 j -1 k=0 β i,j,k ψ j,k (t), t ∈ [0, 1], with J = log 2 (T ), α i,0,0 ≈ √ T 1 0 W i (t)φ 0,0 (t)dt and β i,j,k ≈ √ T 1 0 W i (t)ψ j,k (
t)dt are the empirical wavelet coefficients of the ith individual. A discrete wavelet transform (DWT) corresponds to the computation of these coefficients. In practice, a fast wavelet decomposition and reconstruction algorithm can be computed using the algorithm proposed by [START_REF] Mallat | A theory for multiresolution signal decomposition: the wavelet representation[END_REF] (in only O(T ) operations). As mentioned in the introduction, a simple shift in the observed function will potentially result in a significant change in the DWT. Since we use the latter for feature extraction and the observed curves can start at different times, such behavior is not suitable.

In the TI case, we consider the fast translation-invariant discrete wavelet transform (TIDWT) developed by [START_REF] Coifman | Translation-invariant de-noising[END_REF], in a denoising framework. This transformation has been independently discovered, on several occasions, in different communities, and has received several different names, including the "à trous" algorithm [START_REF] Holschneider | A realtime algorithm for signal analysis with the help of the wavelet transform[END_REF]; [START_REF] Dutilleux | An implementation of the "algorithme à trous" to compute the wavelet transform[END_REF], the undecimated DWT Lang et al. (1996), the shift-invariant DWT Lang et al. (1995) or the stationary DWT Nason and Silverman (1995), to name just a few (see, e.g. [START_REF] Fowler | The redundant discrete wavelet transform and additive noise[END_REF] for a review of some of the various different variants). There are many ways to implement this transformation, and many ways to represent the resulting overcomplete set of wavelet coefficients. We have chosen to focus on the TIDWT of [START_REF] Coifman | Translation-invariant de-noising[END_REF], which provide equivalences and a way to go from one to the other of these representations, for example with the stationary DWT of [START_REF] Nason | The stationary wavelet transform and some statistical applications[END_REF]. The main difference with the orthogonal case is that the dictionary is now a tight frame instead of an orthonormal basis (see (Mallat, 2008, Chapter 5)) and the number of coefficients per scale is no longer dyadic but of length T (see [START_REF] Coifman | Translation-invariant de-noising[END_REF] for more details). This wavelet transform is called translation-invariant by [START_REF] Coifman | Translation-invariant de-noising[END_REF] since the whole dictionary is invariant under circular translation. As in the traditional case, TIDWT is calculated by a series of decimation and filtering operations, only the additional circulant shift S h is added and the corresponding wavelet dictionary is obtained by sampling the locations more finely (i.e., one location per sample point). TIDWT consists of calculating the DWT of the shifted data for each shift h ∈ {0, . . . , T -1}. [START_REF] Coifman | Translation-invariant de-noising[END_REF] propose an algorithm to perform this transformation in O(T log 2 T ) operations (we used the R package rwavelet which provides an implementation [START_REF] Navarro | R package rwavelet: Wavelet analysis[END_REF]). The invariance property of their construction is formally expressed in terms of the circulant matrix containing the wavelet coefficients (see (Coifman and Donoho, 1995, eq. (3))). In other words, for a curve W i translated by h the wavelet coefficients at each scale will be the same up to some permutation. Thus the norm of the latter is preserved scale by scale.

The redundancy of TIDWT makes it possible to detect the presence of hidden information such as stationary or non-stationary patterns as well as their location, making it particularly suitable for clustering purposes. This type of invariant representation has been exploited in many applications (such as denoising [START_REF] Coifman | Translation-invariant de-noising[END_REF] or texture image classification and segmentation [START_REF] Unser | Texture classification and segmentation using wavelet frames[END_REF]). In addition, the use of wavelets allows to compress the information contained in the time series into a small number of wavelet coefficients. Following [START_REF] Antoniadis | Clustering functional data using wavelets[END_REF], we characterize each time series by the vector of the energy contribution of their wavelet coefficients at each scale with the difference that the coefficients are calculated by TIDWT instead of DWT. This extension is possible because the expansion being in a tight frame, the norm is also conserved (see (Mallat, 2008, Chapter 11) for more details). More precisely, using Parseval's identity, we have

W i 2 2 = 2 -J T -1 k=0 α 2 i,0,k + J j=1 2 -j T -1 k=0 β 2 i,j,k = 2 -J θ i0 2 2 + J j=1 2 -j θ ij 2 2 , (4)
where θ ij = (α i,0,0 , . . . , α i,0,T -1 , β i,j,0 , . . . , β i,j,T -1 ) and the factor 2 -j is used to compensate for the redundancy of this representation. The global energy W i 2 2 of W i is decomposed into a small number of components. The representation (4) is made of the components of the invariant version of the discrete wavelet scalogram (as defined in [START_REF] Arino | Wavelet scalograms and their applications in economic time series[END_REF]) can be considered as the TIDWT analogue of the well-known periodogram of the spectral analysis of time series. Similar to how the periodogram produces an ANOVA decomposition of a signal's energy into different Fourier frequencies, the scalogram breaks down the energy into "level components" that indicate at which levels of resolution the energy of the observed function is concentrated. A function that is relatively smooth will have most of its energy concentrated in the low levels j, resulting in a θ ij that is large for small j and small for large j. A function with many high-frequency oscillations will have much of its energy concentrated in the highresolution wavelet coefficients. Therefore, how these energy components are distributed and contribute to the overall energy of a signal is the key fact that we will exploit for clustering. Thus, denoting by y ij the log total squared norm at scale j for the ith individual, we have

y ij = ln θ ij 2 2 , ∀j = 0, . . . , J, i = 1, . . . , n.
(5)

Clustering will therefore be carried out on the basis of the log squared norm of the TI wavelet coefficients at each scale. Thus, this criterion is not sensitive to the origin of the curves, so it seems relevant given the nature of the data motivating this work.

Adjustment on the population risk factors

In this section, we consider the regressions of the features extracted by the wavelet decomposition on the population risk factors. The following lemma shows that the noises of these regressions retain the cluster information given by the daily COVID-19 death curves and permit the information of the population risk factors to be considered in the clustering procedure. Thus, the effects of the risk factors can be estimated from the wavelets coefficients with a nonparametric regression (see ( 6)). Then, the cluster memberships of the region can be assessed from the residuals of the regressions (see ( 7)) since this residuals keep the discriminative information contained in the original data. Note that the same nonparametric function is used for the regression of each feature (i.e., j = 0, . . . , J).

Lemma 1. Let data W 1 , . . . , W n arise from (2) and J + 1 features y i1 , . . . , y iJ are defined by (5) from the wavelet decomposition of W i . Let the centered features be defined by

y ij = y ij -∆ j , ∆ j = E[ln µ(X i )] + L =1 π E[ 1 2 ln v j + ε i j 2 2 ]
, and v j and ε i j are computed from u and ε i , respectively (their formal definition is given in the proof of the lemma presented in Appendix).

y ij = m(x i ) + ξ ij , (6)
with m(x i ) = ln µ(x i ) -E[ln µ(X i )], E[m(X i )] = 0 and E[ξ ij ] = 0,
where the covariates x i and the noise of the regression ξ ij are independent,for j = 0, . . . , J. Then, the noises ξ i follow a mixture model with latent variable Z i defined by the density

g(ξ i ) = L =1 π g (ξ i -λ ), (7)
where λ = (λ 1 , . . . , λ J ) ,

λ j = E[ 1 2 ln v j + ε i j 2 2 ] -∆ j + E[ln µ(X i )]
and g 1 , . . . , g are densities of centered distributions.

We consider the single-index regression defined by m(x i ) := ν(x i γ).

(8)

This semiparametric approach is flexible, and avoids the assumptions of the parametric approaches that can be violated and the curse of dimensionality of the full nonparametric approaches. The parameter of the index γ permits population characteristics having a protective or compounding effect to be detected. Moreover, considering two sets of covariates X i and X i , the difference ν(X i γ) -ν(X i γ) can be interpreted as the logarithm of an odd ratio.

The single-index approach requires a methodology for estimating γ and m, with m being in a function space. A common approach, that avoids a simultaneous search involving an infinite-dimensional parameter, is the profiling [START_REF] Severini | Profile likelihood and conditionally parametric models[END_REF][START_REF] Liang | Estimation and testing for partially linear single-index models[END_REF], which defines ν(x i γ) := ν γ with

ν γ (t) = E[Y ij | X i γ = t],
j ∈ {0, . . . , J} and t ∈ R.

Hence, one expects that, for each x i , the true value of the parameter, denoted by γ realizes the minimum of

γ → J j=0 E[{Y ij -ν γ (x i γ)} 2 | X i = x i ]. (10)
However, even if m γ is well defined for any γ ∈ R d , the vector γ is not identifiable and only its direction could be consistently estimated. Thus, there are two common approaches to restrict γ for identification purposes: either fix one component equal to 1 [START_REF] Ma | Doubly robust and efficient estimators for heteroscedastic partially linear single-index model allowing high-dimensional covariates[END_REF], or set the norm of γ equal to 1 and fix the sign of its first component to be positive [START_REF] Zhu | Empirical likelihood confidence regions in a partially linear single-index model[END_REF]. The estimation of the single-index regressions is performed by considering the empirical counterpart of ( 9) and ( 10)

γ = arg min γ n i=1 J j=1 ŷ ij -νγ (X i γ) 2 , ŷ ij = y ij - 1 n n i=1
y ij , and νγ (u) =

1 nh n i=1 ŷ ij K X i γ-u h 1 nh n i=1 K X i γ-u h ,
where K is a kernel and h a bandwidth. The estimation procedure is implemented in the R package regpro [START_REF] Klemela | regpro: Nonparametric regression[END_REF]. The clustering of the regions is also performed on the residuals ξi = ( ξi0 , . . . , ξiJ ) defined by ξij = ŷ ij -νγ (X i γ).

Nonparametric clustering of the regions

A wide range of literature focuses on models assuming that, conditionally on knowing the particular cluster the subject i came from, its features are independent [START_REF] Chauveau | Semi-parametric estimation for conditional independence multivariate finite mixture models[END_REF]. Thus, we consider that the conditional distribution of the ξi given cluster membership is defined as a product of univariate densities. Note that this assumption is standard when nonparametric modelling of the mixture components is wanted. Moreover, this assumption is suitable when J is quite large with respect to n because it limits the number of parameters to be considered. Note that this assumption imposes non-explicit constraints on the distribution of the noises ε i defined in (1). Note that this assumption could be relaxed by in a parametric framework [START_REF] Bouveyron | The discriminative functional mixture model for a comparative analysis of bike sharing systems[END_REF] but also in a nonparametric framework [START_REF] Zhu | Clustering via finite nonparametric ica mixture models[END_REF]. In numerical experiments (see Section 4.1), we illustrate the relevance of the conditional independence assumption. Therefore, the clustering of the region is performed by considering the marginal density defined by

g( ξi ; λ) = L =1 π J j=1 g j ( ξij ), ( 11
)
where λ groups the mixing proportions π 1 , . . . , π L (where π > 0 and L =1 π = 1) and the univariate densities g j . The model ( 11) is identifiable, up to a swapping of the cluster labelling, if the densities g j are linearly independent (see Theorem 8 of [START_REF] Allman | Identifiability of parameters in latent structure models with many observed variables[END_REF]). Considering a multivariate kernel defined as a product of J univariate kernels K, the maximum smoothed log-likelihood estimator λ (MSLE) is obtained by maximizing the smoothed log-likelihood (λ) [START_REF] Levine | Maximum smoothed likelihood for multivariate mixtures[END_REF], such that λL = arg max λ (λ; L) and

(λ; L) = n i=1 ln L =1 π J j=1 N g j ( ξij ) ,
where

N g j ( ξij ) = exp Ω j 1 h j K ξij -u h j ln g j (u)du ,
and h 1 , . . . , h J are the bandwidths (i.e., h j > 0 and h j = o(1) for j = 1, . . . , J). Considering the MSLE is more convenient than considering the maximum likelihood estimate because the MSLE can be obtained by a Majorization-Minimization algorithm (see [START_REF] Levine | Maximum smoothed likelihood for multivariate mixtures[END_REF] for details on the algorithm and [START_REF] Zhu | Theoretical grounding for estimation in conditional independence multivariate finite mixture models[END_REF] for recent developments) implemented in the R package mixtools [START_REF] Benaglia | mixtools: An r package for analyzing finite mixture models[END_REF].

Clustering is achieved by computing the MSLE because this estimator permits a soft assignment where the conditional probability that subject i belongs to cluster , denoted by t i ( λ), can be obtained

t i ( λL ) = π J j=1 N ĝ j ( ξij ) L =1 π J j=1 N ĝ j ( ξij )
.

Moreover, a hard assignment can be achieved by applying the maximum a posteriori rule (leading that ẑi = 1 if = arg max t i ( λ) and ẑi = 0 otherwise).

4 Numerical experiments

Investigating the strengths of the proposed approach

Data are independently generated from (1) with T = 512, K = 3, unequal proportions (π 1 = 0.5, π 2 = 0.25 and π 3 = 0.25) by considering two scenarios:

• For scenario 1, we define

u (δ i ) (t) = r (t -δ i )1 {r (t-δ i )>0} if t > δ i 0 otherwise ,
with r 1 (t) = sin(2.5πt/T ), r 2 (t) = 1.5 sin(2.5πt/T ) and r 3 (t) = sin(2πt/T ) and the vector of noises ε i is composed of independent Gaussian random variables with mean 0 and variance 1/2.

• For scenario 2, we define

u (δ i ) (t) = r (t -δ i ) if t > δ i 0 otherwise ,
and consider the same simulation setting presented in (Antoniadis et al., 2001, Section 4.1). That is, r 1 (t) = sin(2πt/T )+sin(5πt/T ) and for the second and third clusters, we use two strictly stationary autoregressive Hilbertian (ARH) processes with the same parameters (see (Antoniadis et al., 2001, Section 4.1) for more details).

For both scenarios, the time shift δ i of curve i is generated from a uniform distribution defined on {0, . . . , ∆}. Moreover, we generate two covariates from standard Gaussian distributions independently, and we consider γ

= (1/ √ 2, 1/ √ 2) and ν(t) = (1 + qt 2 )/(1 + q).
Relevance of the translation-invariant feature extraction This experiment illustrates the importance of using translation-invariant decomposition when the observed curves are time-shifted. Thus, data are generated with different values of ∆ and no covariate effects (i.e., q = 0). To investigate the relevance of the proposed approach, we compare the translation-invariant wavelet decomposition (TI-wave) with different functional basis: polynomial basis (Poly-basis), Bsplines (Bsplines) and exponential basis (Exp-basis). Each basis decomposition considers 9 basis elements. For each basis decomposition, the clustering is achieved by fitting the nonparametric mixture model defined by ( 11) on the basis coefficients. Figure 2 shows the Adjusted Rand Index (ARI) obtained on 100 replicates by the different basis decomposition for different sample sizes and different values of ∆, under both scenarios. When there is no time shifting (i.e., ∆ = 0), under scenario 1 the polynomial basis outperforms all the other approaches, while the translation-invariant wavelet basis, polynomial basis and the exponential basis perform similarly under scenario 2. Note that exponential basis with 9 elements seems not to be relevant for these scenarios. When the time-shifting increases, then all the results obtain by the different basis are (strongly) deteriorated expect those obtained by the translation-invariant wavelet frame. These results illustrate the importance of applying a translation-invariant clustering method to analyze the COVID-19 data. In the next section, we compare the proposed approach with other translation-invariant approaches.

Relevance of the covariate adjustment This experiment illustrates the importance of considering the covariate effects. Thus, data are generated with different covariate effects (q ∈ {0, 0.1, 0.2, 0.3}) and no time shifting (i.e., ∆ = 0). To investigate the relevance of the proposed approach, we compare the results obtained by the proposed three-step approach with those obtained by neglecting the covariate effect (i.e., ν is supposed to be the constant function equals to one). Figure 3 presents the distribution of the ARI obtained by the proposed approach that considers the covariate effect and by a clustering performed on the translation-invariant wavelet frame which neglects the covariate effect. When there is no covariate effect (i.e., q = 0), the approach modeling the covariate effects obtained results that are slightly outperformed by the approach neglecting the covariate effects. However, when there is a covariate effect, the results obtained by the approach neglecting the covariate effects are strongly deteriorated while modeling the covariate effect allows the partition to be consistently estimated.

Figure 3: Boxplots of the Adjusted Rand Index obtained on 100 replicates by considering (proposed) and by neglecting (No cov) the covariate effects when data are generated with different covariate effects q (0, 0.1, 0.2 and 0.3).

Relevance of the nonparametric clustering This experiment focuses on the influence of the mixture model used for clustering. Thus, data are generated with a time shift defined by ∆ = 20 and a covariate effect defined by q = 0.1. To investigate the relevance of the mixture model defined by ( 11), clustering of ξ1 , . . . , ξn is performed by considering (11), the semiparametric mixture considering the intra-component dependencies proposed by [START_REF] Zhu | Clustering via finite nonparametric ica mixture models[END_REF], and the parsimonious Gaussian mixture models implemented in the R package HDclassif [START_REF] Bergé | Hdclassif: An r package for model-based clustering and discriminant analysis of high-dimensional data[END_REF]. Note that performing clustering on the basis coefficients with the Gaussian mixture implemented in R package HDclassif is the method proposed by [START_REF] Bouveyron | The discriminative functional mixture model for a comparative analysis of bike sharing systems[END_REF] and implemented in the R package funFEM [START_REF] Bouveyron | The discriminative functional mixture model for a comparative analysis of bike sharing systems[END_REF]. Figure 4 presents the distribution of the ARI according to the clustering method. Overall, the results of the three clustering methods are similar, despite the fact that the results obtained by considering intra-component dependencies slightly deteriorate the results of the nonparametric mixture. Under scenario 1, the proposed method (nonparametric mixture model with components defined as a product of univariate densities) slightly outperforms the parsimonious Gaussian mixture. Under scenario 2, all the methods perform similarly. These results show that the choice of the clustering method (step 3 of the proposed approach) is less impacting that the choice of the basis expansion (step 1 of the proposed approach) and the modeling of the covariate effect (step 2 of the proposed approach).

Figure 4: Boxplots of the Adjusted Rand Index obtained by the different mixture models (Sparse non-parametric corresponds to the proposed approach, HDclassif and a Full nonparametric mixture) on 100 replicates according to the sample size.

Comparing the proposed approach with other methods used for COVID-19 studies

In this experiment, we compare the proposed three-step approach with alternatives used in other studies about COVID-19. Thus, data are generated from both scenarios described in the previous section with different values of ∆ and q. Figure 5 presents the distribution of the ARI obtained by the proposed method, by the non-negative matrix factorization (NMF) approach [START_REF] Chen | Clustering us states by time series of covid-19 new case counts with non-negative matrix factorization[END_REF] and by using the similarity matrix given defined by the DTW approach for a clustering achieved via spectral clustering [START_REF] Allem | A spectral clustering approach for the evolution of the covid-19 pandemic in the state of rio grande do sul, brazil[END_REF] and via hierarchical ascendant classification [START_REF] Park | Clustering county-wise covid-19 dynamics in north carolina[END_REF] defined with different criteria (complete linkage, ward criterion, single linkage). When there is no time shifting and no covariate effect (∆ = 0 and q = 0), NMF outperforms the proposed method under both scenarios and the spectral clustering performed on the DTW dissimilarity matrix outperforms the proposed method under scenario 1. However, our real data analysis considers regions from America and Europe, then time-shifting occurs. The results of the NMF methods are strongly deteriorated when data are time-shifted (∆ = 60) under both scenarios. Thus, NMF seems to be relevant to focus on a particular region where the diseases appear in the same time. Moreover, when there is a covariate effect, the results obtained by DTWbased approach deteriorated. To the best of our knowledge, there is no method that allows covariate to be included in DTW-based similarity matrix. Thus, the proposed method seems to be the single approach that allows to deal with time-shifting and to consider covariate adjustments.

Figure 5: Boxplots of the Adjusted Rand Index obtained on simulated data by the different clustering methods used in COVID-19 studies.

We now use the COVID-19 data to generate the samples. We consider the daily COVID death curves of Germany, Denmark and France (denoted by W1 , W2 and W3 ). Data are independently generated from (1) with T = 512, K = 3, equal proportions and

W i (t) = L =1 z i µ(x i γ)( W (δ i ) + ε (δ i ) i ).
The time shift δ i of curve i is generated from a uniform distribution defined on {0, . . . , ∆}, the two covariates are generated from standard Gaussian distributions independently, γ = (1/ √ 2, 1/ √ 2) , µ(t) = (1 + qt 2 )/(1 + q) and the noises ε i arise independently from a centered Gaussian distribution with variance v 2 . Figure 6 presents the distribution of the ARI obtained by the competing methods. It shows that the proposed method outperforms the other competing methods when there is time shifting and a covariate effect.

Investigating the robustness of the proposed approach

We now investigate the robustness of the proposed approach to spatial dependencies. Indeed, these dependencies are not considered during the estimation but can occur in the COVID-19 data set. Thus, we consider the two experiments described in Section 4.2 but relaxing the assumption of independence between the ε i . Thus, we now consider that the each vector (ε 1 (t), . . . , ε n (t)) independently arises from a centered Gaussian random field with parameters σ 2 and φ where the distance between two sites is computed from their spatial coordinates that are generated uniformly in [0, 1] 2 . Figure 5 presents the distribution of the ARI obtained by the competing methods when the noises follow Gaussian random fields. Figure 8 presents the distribution of the ARI obtained by the competing methods when the data are simulated from the COVID-19 dataset. Results show that the proposed method outperforms the other competing methods when there is time shifting and a covariate effect. Moreover, the spatial dependency does not deteriorate the results. of the curve due to extreme points). Overall, the shape of μ implies that the more the index is, the more is the COVID-19 death rate. Moreover, the estimator of the density of the index is smoothed enough to ensure a consistent estimation of the regression function. Note that the three relevant covariates in X i are composed of one PCA-scores measuring environmental risk factors and two PCA-scores measuring the medical risk factors. Thus, these results confirm that diabetes, overweight, smoking, pulmonary disease, cardiovascular disease, HIV, hypertension and age are risk factors for the COVID-19 (see Table 2 to interpret the relation between the PCA-scores and the risk factors). These results are in agreement with the main mortality risk factors identified in medical publications (see Zhou et al. (2020, Table 1) and Gupta et al. (2020, Figure 2)). Moreover, the population density and the concentrations of nitrogen dioxide and fine particulate matter also increase the COVID-19 mortality risk. Additionally, these results align with the findings of other works [START_REF] Sy | Population density and basic reproductive number of covid-19 across united states counties[END_REF][START_REF] Copat | The role of air pollution (pm and no2) in covid-19 spread and lethality: a systematic review[END_REF][START_REF] Pozzer | Regional and global contributions of air pollution to risk of death from covid-19[END_REF].

Clustering of the regions

Clustering is performed by considering a number of clusters between one and ten. Selecting the number of clusters is still a difficult task in nonparametric mixtures. Indeed, despite recent works [START_REF] Kasahara | Non-parametric identification and estimation of the number of components in multivariate mixtures[END_REF][START_REF] Kwon | Estimation of the number of components of non-parametric multivariate finite mixture models[END_REF] presenting elegant methods for selecting the number of clusters based on the constraint of linear independence between the univariate densities required for model identifiability [START_REF] Allman | Identifiability of parameters in latent structure models with many observed variables[END_REF], these Table 4: Mean of the DTW distances between the adjusted curves of two classes.

• Canada, Germany, Denmark, Dominican Republic, Estonia, Greece, Haiti, Netherlands, Norway, Serbia, Florida, Indiana and Ohio.

• Cuba, Cyprus, Guyana, Iceland, Jamaica, Latvia, El Salvador, Alaska, Delaware, District of Columbia, Hawaii, Idaho, Illinois, Louisiana, Maine, Michigan, Minnesota, Mississippi, Missouri, Nebraska, Nevada, New Hampshire, New Mexico, North Carolina, Oregon, Pennsylvania, Tennessee, Texas, Utah, Vermont, Washington and Wisconsin.

• Chile, Spain, Ireland, Sweden, Alabama, Arkansas, Colorado, Georgia, Iowa, Kansas, Kentucky, Maryland, Montana, North Dakota, Oklahoma, South Carolina, South Dakota, Virginia, West Virginia and Wyoming.

• Austria, Belgium, Bulgaria, Brazil, Switzerland, Colombia, Costa Rica, Czechia, France, United Kingdom, Croatia, Hungary, Italy, Panama, Peru, Poland, Portugal, Paraguay, Romania, Slovenia, Slovakia, Arizona, California, Connecticut, Massachusetts, New Jersey, New York, Rhode Island and Uruguay.

We now describe the clusters based on summary statistics presented in Table 5 and the curves adjusted with the population risk factors (W i /μ(X i γ)). Figure 11 presents the curves normalized by the population risk effects W i /μ(X i γ) per cluster to illustrate our cluster interpretation. In addition, Figure 12 presents the original curves per clusters to illustrate the impact of the normalization by the population risk effects. The labels of clusters are ordered by the impact of the COVID-19 during the studied period. Thus, we notice that cluster 1 has the smallest mean of COVID-19 daily death rates over the studied period and cluster 5 has the highest mean. These results were expected because the distribution of μ(X i γ) is supposed to be the same among clusters. Clusters can be interpreted, as follows:

• Cluster 1 contains 13 regions that are mainly unaffected from the disease.

• Cluster 2 contains 32 regions that suffer from multiple small waves of deaths.

• Cluster 3 contains 20 regions mainly affected by one major wave with steep increasing and decreasing rates.

• Cluster 4 contains 29 regions strongly impacted by the disease and that suffered for two waves at the minimum. In addition, for the regions that suffered for at least three waves of death, there is a wave of deaths which starts before the previous wave ends.

Clusters analysis example: disparities and policy decisions

Amid unforeseen difficulties related to COVID-19, policymakers have resorted to various interventions in attempts to curb the spread of the corona virus. As seen in Sections 5.1 and 5.2, by adjusting the risk factors, the proposed approach allows us to detect vulnerable regions and compare them fairly. Additionally, the homogeneity within cluster enables us to analyze the disparities in factors between clusters; for instance, in this paper, we are looking at the government responses. Note that one may also be interested to study COVID-19's impact on the population mental health or the economic indicators. Figure 13(a) shows the distribution of the value for each indicator, i.e., containment health, government response and stringency, with respect to its clusters. We can observe that the indicators take larger values for regions of assigned in cluster 4 (containing the most impacted regions). The same phenomenon is observed for the specific measures adopted by the governments presented in Figure 13(b).

Cluster 4, which contains the regions strongly affected by the virus, seem to take more stringent measures than the other regions. One plausible explanation is that the policymakers may abruptly enforce restrictive countermeasures in hopes to lower the increasing mortality rate, thus the observed rapid descent. Hence, the effectiveness of government response may depend on the timing of the measure's implementation, the duration and the stringency [START_REF] Cheng | Covid-19 government response event dataset (coronanet v.1.0)[END_REF][START_REF] Haug | Ranking the effectiveness of worldwide covid-19 government interventions[END_REF]. To further understand the relationship between the mortality rate and the various policies adopted in reaction to the COVID-19, a more precise analysis can be achieved by fitting a model for multivariate non-stationary time series per clusters (the time series being the daily death rate and the index mentioned above). Thus, a description of the relation between the government responses and the COVID-19 death curve can be done using the model parameters. This can be achieved by modelling the dependencies between the time series [START_REF] Molenaar | Dynamic factor analysis of nonstationary multivariate time series[END_REF][START_REF] Sanderson | Estimating linear dependence between nonstationary time series using the locally stationary wavelet model[END_REF] or by a multiple change-point detection [START_REF] Cho | Multiple-change-point detection for high dimensional time series via sparsified binary segmentation[END_REF]. Moreover, a comparison of the government interventions to COVID-19 can be conducted by comparing the models fitted for each cluster.

Conclusion

This paper introduces CRFTIW, a new method to cluster functional data when observations are shifted and external covariates are allowed to have a scale effect. CRFTIW is a three-step approach developed for the purpose of analyzing geographical disparities of the COVID-19 impact (measured by the daily number of deaths per million people).

As a first step of CRFTIW, feature extraction is performed with TI wavelets. While providing an adapted and compact representation of the data, it also permits us to deal with the different times of arrivals of the disease. The main limitation of this approach lies in the dyadic data constraint of the considered sample. This issue could be overcome, for example, by using second generation wavelets and in particular the lifting scheme [START_REF] Sweldens | The lifting scheme: A construction of second generation wavelets[END_REF]), but would lose the property of translation-invariance. However, extending this construction, while preserving both the property of translation-invariance and the property of conservation of the norm of the coefficients, seems to be an open question that we leave for future work.

As a second step of CRFTIW, the effect of the population risk factors on the extracted feature is estimated and thus regions can be compared as if they have the same sensitivity of the population to the disease. This step is crucial because we aim to investigate the impact of policy decisions. Obviously, if the purpose is to investigate geographical disparities of the impact of the disease, then no adjustment on the population risk factors should be considered. In such case, CRFTIW can still be used by considering µ(X i γ) = 1. In the analysis of COVID-19, we consider four main factors of comorbidity (overweight, diabetes, chronic pulmonary and kidney diseases). These factors are well-known to increase the risk of COVID-19 mortality, [START_REF] Zhou | Clinical course and risk factors for mortality of adult inpatients with covid-19 in wuhan, china: A retrospective cohort study[END_REF][START_REF] Gupta | Factors associated with death in critically ill patients with coronavirus disease 2019 in the us[END_REF] and Wilks' theorem for the empirical likelihood permits us to conclude on the significance of the estimated coefficients. However, we advise considering factors already known to be compounded for the disease, and not to use CRFTIW to investigate the impact of population risk factors. Indeed, CRFTIW does not perform variable selection of the population risk factors and does not permit concluding on causality of the factors.

As a third step of CRFTIW, a nonparametric mixture is used to achieve the clustering with the assumption that the density of the component is defined as a product of univariate densities. Numerical experiments presented in the paper suggested that considering a more complex model deteriorated the results when the sample size is small (like in the COVID-19 application). However, if the data to be analyzed are composed of several observations, more advanced models could be used [START_REF] Mazo | Constraining kernel estimators in semiparametric copula mixture models[END_REF][START_REF] Zhu | Clustering via finite nonparametric ica mixture models[END_REF].

Through the COVID-19 dataset, we had illustrated the importance of adjusting the population risk factors, allowing us to compare regions with a 'standard' comorbidity. Thus, CRFTIW found five clusters justified by the mortality rate and curvature. Regions within clusters are varied geographically with different onsets, validating the property of translationinvariance of the proposed method. In addition, as we illustrated, investigations on the effectiveness and agility of government response, the consequences on economic indicators or the impact on human mental health, could be achieved by studying disparities of the indicators between clusters.

Despite the model being translation invariant, the time between the arrivals of two waves is discriminative. We argue that this time is important; it determines whether the health facilities have any breaks between waves. Indeed, countries suffering for successive waves of COVID-19 have to postpone non-emergency surgical operation or early cancer detection. Thus, using the proposed clustering, we could investigate the impact of COVID-19 on the global quality of care. Note that an alternative clustering approach could focus only on the death peaks thus neglecting the time between waves. In such a case, the proposed approach is not suitable and we advice to use time scaled clustering [START_REF] Tang | Time-synchronized clustering of gene expression trajectories[END_REF].

Driven by the COVID-19 dataset, we developed this novel approach. However, its application is not limited only to COVID-19. For instance, the problem of time-shifts is also observed in electrocardiogram heartbeat, which [START_REF] Annam | Time series clustering and analysis of ecg heart-beats using dynamic time warping[END_REF] tackles when clustering heartbeat abnormalities. Our approach could not only handle the time-shift issue, but also allow for adjusting with regards to plausible factors that may influence the heartbeat. Further, this could extend beyond medical settings: in motion capture [START_REF] Li | Time series clustering: Complex is simpler! In ICML[END_REF], there is interest in categorizing types of motion. This could be used for fitness applications to identify whether a person is running or walking, where the motions may begin at different times on separate observed sequences showing a need for the TI property. Since motion can come from different participants, covariate adjustment could also be beneficial for such data.

Figure 1 :

 1 Figure 1: Illustration of different arrivals of the COVID-19 among the different regions (Europe, North America, Latina America).

Figure 2 :

 2 Figure 2: Boxplots of the Adjusted Rand Index obtained on 100 replicates according to the basis decomposition (TI-wave, Poly-basis, Bspline and Exp-basis), the sample size (50, 100 and 200) and the value of ∆ (0, 20, 40 and 60).

Figure 6 :

 6 Figure 6: Boxplots of the Adjusted Rand Index obtained on data simulated from the COVID-19 dataset by the different clustering methods used in COVID-19 studies.

Figure 7 :

 7 Figure 7: Boxplots of the Adjusted Rand Index obtained on simulated data by the different clustering methods used in COVID-19 studies when the noises follow Gaussian random fields.

Figure 8 :

 8 Figure 8: Boxplots of the Adjusted Rand Index obtained on data simulated from the COVID-19 dataset by the different clustering methods used in COVID-19 studies when the noises follow Gaussian random fields.

Figure 9 :

 9 Figure 9: Function μ modeling the effect of the index on the COVID-19 death curves (on the left) and density of the index X i γ{2,5} (on the right).

Figure 10 :

 10 Figure 10: First PCA map computed from the ξ 1 , . . . , ξ n . Colors indicate the estimated partition in five classes (red: class1; maroon: class2; green: class3; blue: class4; purple: class5.

Figure 11 :

 11 Figure 11: Curves W i /μ(X i γ) displayed per cluster (1 to 4 from the left to the right).

Figure 12 :

 12 Figure 12: Original curves W i displayed per cluster (1 to 4 from the left to the right).

Figure 13 :

 13 Figure 13: Boxplot of the overall indicator values for each cluster (a) and for specific measures adopted for each cluster (b): School closing (C1), Workplace closing (C2), Cancel public events (C3), Testing policy (H2) and Contact tracing (H3).

Table 5 :

 5 Statistics per cluster cluster proportion Normalized deaths UnNormalized deaths Covariate effect

		π	mean	sd.	mean	sd.	mean	sd.
	Class 1	0.14	860.10	397.02	958.87	638.69	1.01	0.41
	Class 2	0.33	1030.55	414.88	1219.62	711.38	1.11	0.36
	Class 3	0.22	1463.05	364.75	1736.15	350.80	1.23	0.28
	Class 4	0.31	2532.97	961.54	2253.94	888.82	0.94	0.32

Acknowledgements

The authors thank the editor, the associate editor, and two referees for their many constructive comments and suggestions, which greatly improved the quality of the article.

Investigating geographical disparities for COVID-19

The proposed approach allows risk factors to be taken into account, and thus to investigate the geographical disparities of the COVID-19 that are not explained by geographical disparities of the population risk factors. Thus, this adjustment permits similarities and disparities of the disease impacts to be detected, conditional on the population characteristics. In Section 5.1, we interpret the impact of the population risk factors on the COVID-19 death curves and we confirm this interpretation with literature review. In Section 5.2, we describe the clustering outputs. Finally, Section 5.3 presents an a posteriori comparison of the clusters with respect to different policy decisions (e.g., lockdown characteristics).

Population risk factors

During the application, the first principal component obtained on the PCA of the environmental variables and the first two principal components obtained on the PCA of the medical variables can be considered. Indeed, the second principal component obtained on the PCA of the environmental variables and the third principal component obtained on the PCA of the medical variables are not relevant in the single-index regression (see ( 8)). The relevance of the estimated coefficients of the single-index regression γ is investigated by using empirical likelihood (Du Roy de Chaumaray et al., 2021). Table 1 presents the estimator of the regression coefficients and its p-value provided by using the empirical likelihood ratio. It also presents the estimators of the coefficients obtained by assuming that one of the four variables does not have an effect on the target variables (where

under the constraint γ {ω} j = 0 for any j ∈ ω where X i ∈ R 5 is the vector composed of the first two principal component obtained on the PCA of the environmental variables and the first three principal component obtained on the PCA of the medical variables).

Risk factor

Estimators γ {2,5} γ {2,4,5} γ {2,3,5} γ {1,2,5} The impact of these covariates on the daily death rate is made throughout the index X i γ{2,5} where γ{2,5} = (0.533, 0.000, 0.658, 0.531, 0.000) and the regression function μ. Figure 9 shows the estimator of the μ and the density of the index for a range covering 90% of the observed index (this trimming, performed only for this plot, avoids over-interpretation ] and the bounds of the bins are defined by the empirical quantiles of ξij . The approach permits to select K = 4 clusters thus corresponding in the first an elbow in the values of the smoothed log-likelihood that are presented in Table 3. Note that a second elbow occurs and suggests that a partition into seven clusters is also of interest. In the following, we provide the interpretation of the clustering output into 4 clusters. Figure 10 represents the regions in the first PCA map computed from the ξ 1 , . . . , ξ n and it indicates the partition with colors (Table 6 presented in Appendix B gives the relation between the ID of the regions and their names). Hence, it shows that the estimated partition provides no overlapping cluster in the space of the translation-invariant wavelet frame adjusted with the population risk factors. Moreover, to investigate the within-class homogeneity between the curves adjusted by the population risk factors, Table 4 presents, for each class k and , the mean of DTW distances between each curve belonging to class k and each curve belonging to class (i.e., {i:

A Proof of Lemma 1

Proof of Lemma 1. We consider the matrix representation for the discrete wavelet transform. For the orthogonal case, the discrete transform of a vector w ∈ R T is represented by an T ×T orthogonal matrix M (see, e.g., [START_REF] Mallat | A wavelet tour of signal processing: the sparse way[END_REF]).

In the translation-invariant case, the downsampling of the locations k is discarded at each scale j, and the decomposition can thus be written as follows

The matrix representation is therefore given by an (J +1)T ×T matrix M which corresponds to the row-wise concatenation of the J, T × T matrices which yield the wavelet coefficients at scale j plus the matrix which produces the scaling coefficients (see, e.g., [START_REF] Berkner | Smoothness estimates for soft-threshold denoising via translation-invariant wavelet transforms[END_REF]; [START_REF] Coifman | Translation-invariant de-noising[END_REF] for a full representation, in terms of filters). The representation at scale j for W i given z i is defined by M j W i where M j denotes either the T × 2 j or the T × T wavelet transform matrix at scale j, for j = 0, . . . J. Thus, we have

Therefore, the logarithm of the norm at scale j for W i given z i is given by

where

Considering the centered version of the elements leads to (6). Noting that the noise of the regression follows a mixture model where the latent variable is z i (i.e., the same latent variable that the one used to define the marginal distribution of W i ), permits to conclude that clustering the noise of the regression obtained for the different scale is relevant to cluster the curves W i . Finally, noting that by assumption X i is independent to Z i and ε i and that ε i j is defined as a combination of the original ε i permits to conclude on the independence between X i and ε i j .

B Association between the region IDs and their names