Smooth copula-based generalized extreme value model and spatial interpolation for extreme rainfall in Central Eastern Canada

Received: 24 November 2021 Revised: 30 November 2022 Accepted: 18 January 2023 DOI: 10.1002/env.2795 RESEARCH ARTICLE Smooth copula-based generalized extreme value model and spatial interpolation for extreme rainfall in Central Eastern Canada Fatima Palacios-Rodriguez1Elena Di Bernardino2Melina Mailhot3 1Departamento de Estadística e Investigación Operativa, Facultad de Matemáticas, Universidad de Sevilla, Seville, Spain 2Laboratoire J. A. Dieudonné, UMR CNRS 7351, Université Côte d’Azur, Nice, France 3Department of Mathematics and Statistics, Concordia University, Montréal, Quebec, Canada Correspondence Fatima Palacios-Rodriguez, Departamento de Estadística e Investigación Operativa, Facultad de Matemáticas, Universidad de Sevilla, Calle Tarfia sin número, 41012 Seville, Spain. Email: [email protected] Funding information Agence Nationale de la Recherche; Mitacs Globalink; Simons Foundation Abstract This paper proposes a smooth copula-based Generalized Extreme Value (GEV) model to map and predict extreme rainfall in Central Eastern Canada. The considereddatacontainsalargeportionofmissingvalues,andoneobservesseveralnonconcomitant recordperiodsat differentstations.Theproposedtwo-step approach combines GEV parameters’ smooth functions in space through the use of spatial covariates and a flexible hierarchical copula-based model to take into account dependence between the recording stations. The hierarchical copula structure is detected via a clustering algorithm implemented with an adapted version of the copula-based dissimilarity measure recently introduced in the literature. Finally, we compare the classical GEV parameter interpolation approaches with the proposed smooth copula-based GEV modeling approach. KEYWORDS Copula-based clustering, extreme value theory, hydrology, missing values, nonconcomitant record periods, spatial interpolation 1INTRODUCTION Heavy rainfall can have disastrous consequences on health of financial systems and well-being of communities, buildings, infrastructures, transportation systems, and public safety. In practice, the interest is, among others, from a national safety, risk management, and insurance perspectives. For researchers, as Extreme Value Theory (EVT, see, e.g., Beirlant etal., 2004;de Haan& Ferreira,2006) isan areawithgreatrecentinnovativeresults,precipitationlevelsareveryinteresting.Throughouttheyears,floodeventsbecameimportant,bothfrompracticalandtheoreticalperspectives.Forexample, in Canada, it is recent that insurance companies offer flood protections. Before 2013, homeowners would rely on the Disaster Financial Assistance program offered by the federal, provincial, and territorial governments. Now, insurance products are available, and several resources are dedicated to improving flood mapping and mitigation efforts. Statistical models are now adapted with today’s knowledge on extreme events, in order to assess risks depending on precipitation appropriately.Inotherwords,itisdesirabletosetasidesafetycapitalaccordingtoarandomvariableevaluatedasprecisely as possible in terms of distribution, measurements and variability. Extensive literature exists on the spatial mapping or spatial interpolation of extreme rainfall and the approaches are essentially divided in two groups. The first one is mainly composed of spacial interpolation techniques to estimate This is an open access article under the terms of the Creative Commons Attribution-NonCommercial-NoDerivs License, which permits use and distribution in any medium, provided the original work is properly cited, the use is non-commercial and no modifications or adaptations are made. © 2023 The Authors. Environmetrics published by John Wiley & Sons Ltd. Environmetrics. 2023;34:e2795. wileyonlinelibrary.com/journal/env 1of30 https://doi.org/10.1002/env.2795 2of30 PALACIOS-RODRIGUEZ et al. precipitation quantities from marginal extreme value distributions in order to provide return level maps. This classical approach is frequently used, for example, in Beguería and Vicente-Serrano (2006), Cooley et al. (2007), and Kohnová et al. (2009). A comparison of different traditional interpolation methods (inverse distance weighting, nearest neighbor, and kriging) for mapping extreme precipitation in central Slovakia can be found in Szolgay et al. (2009). Similar study is performed in Blanchet and Lehning (2010), or in Das et al. (2020), for mapping snow depth return levels. In Hwang et al. (2012), a two step procedure is proposed where a regression and hydrological models are used to interpolate precipitation. Pointwise return levels, in the Cévennes–Vivarais region (southern part of France), based on nearest neighbor estimators, are obtained in Gardes and Girard (2010). The second group in the extreme return level maps literature is based on the direct estimation of the spatial extremal distribution. Spatial extreme distributions have received a lot of attention in recent years and they represent a well-founded approach which are theoretically preferred to any interpolation method. Several techniques have been developed for the direct estimation of spatial extreme distributions, which involve, among others, extreme-value copulas (see e.g., Joe, 1994; Saad et al., 2015), max-stable processes (see e.g., Davison et al., 2012; Huser et al., 2019;Padoan etal., 2010; Ribatetet al., 2012;Reich& Shaby, 2012; Schlather, 2002; Smith,1990) and Bayesianhierarchicalmodels (see e.g., Johannesson et al., 2021). Notice that the property of max-stablity is classically included in several papers on spatial interpolationofextremes(seee.g.,Blanchet&Davison,2011,forextremesnowdepthmodels).Anextremerainfallmodel assuming regional dependence is proposed in Yoon et al. (2015) and random fields approaches are presented in Sang and Gelfand (2009,2010). Thisarticlefocusesonconstructingthespatialmappingofmaximumprecipitation,usingtheEVTframeworkfor24-h duration rainfall annual maxima in Central Eastern Canada. Despite the fact that several authors have brought efforts in order to provide spatial extreme models for precipitation, the considered dataset used in this article presents at least two interesting aspects which need to be addressed carefully. Inthis paper,wefirst focuson modelingextremerainfallfor 116recordingstationslocatedintheprovinceof Quebec, Nova Scotia, New Foundland, New Brunswick and Prince Edward Island. This is a vast region, with a small proportion (116 stations for more than 1,000,000 km2) of recording stations compared to Lehmann et al. (2016), for example, where 872 and 1348 stations are available on a 156,000 and 580,000 km2regions, respectively, and the literature is quite scarce for this specific area. Note that the size of the meteorological stations network in Canada is a major well-known issue already raised by the Canadian Standards Association. Khedhaouiria et al. (2020) consider another Canadian dataset which focusses on a southern Canadian region. Even given the scarce aspect of the dataset, the obtained results show a robust performance with the considered smooth Generalized Extreme Value (GEV) distribution fitting methods. The same dataset has been analyzed in Perreault et al. (2022) where a Bayesian hierarchical interpolation model is proposed with the spatial effect modeled via Gaussian Markov random fields (see, e.g., Rue & Held, 2005). Results of Section 4can be compared with those of Perreault et al. (2019). Also, we consider a dataset which presents several missing values and nonconcomitant record periods at different stations. The interested reader is referred to Figure E1 (see Data S1) for a graphical illustration of this crucial point. It impliesthatwhenoneaimstomodelthejointbehavioroftheextremerainfallataspecificstation,theestimatedmarginal distribution uses the complete series for a given station, but the copula function representing the dependence (see, e.g., Nelsen, 1999) is only based on the time period where all series were recorded simultaneously. Here, we are facing the problem of estimating parametric multivariate models when unequal amounts of data are available on each variable (see e.g., Patton, 2006). To cope with this situation, we consider the hybrid copula model (see e.g., Segers, 2015). This model is an extension of the classical empirical copula, obtained by combining an estimator of a multivariate cumulative distribution with estimators of the marginal cumulative distributions. Note that, in the missing data framework, the considered marginal estimators are not necessarily equal to the margins obtained through the joint estimator. The main contribution of the present work is the proposed spatial smooth GEV model, mixing response surfaces for the GEV parameters’ models, with a flexible joint dependence framework via a hierarchical copula, taking into account the spatiotemporal dependence structure between the recording stations. The proposed model for the return level maps estimation is based on several inference steps, in order to deal with the challenges of the data, concisely described in the following.Theinterestedreaderis referredtoSection2forthedatadescriptionandSection3.4forthe associateddetailed algorithms. First, in order to handle the dependence between stations, we estimate our adapted version of the copula-based dissimilarity measure recently proposed in Disegna et al. (2017) (see Equation 14). For each couple of stations, this measure is based on a convex combination between the spatial distance and the copula behavior. Pseudo-observations for the 1099095x, 2023, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/env.2795 by Universidad De Sevilla, Wiley Online Library on [19/09/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License PALACIOS-RODRIGUEZ et al. 3of30 copula are built by using the point-wise GEV quantiles with parameters estimated via the L-moment method. Then, on thesepseudo-observationsseriesofdifferentlengthsandnonconcomitantyears,weconsistentlyestimatethedissimilarity by using the hybrid empirical copula (see Algorithm 1). Second, we build a hierarchical copula model via a PAM clustering algorithm based on the estimated dissimilarity, considering the spatial dependence between stations (see Algorithm 2). Finally we include the GEV smooth surface parameters’ with polynomial regression and spline models and we maximize the obtained hierarchical loglikelihood (see Equations 9and 10). Considered covariates included in the considered polynomial/spline models are presented in Table 1. The resulting smooth functions for the GEV parameters are used to provide the estimated return level maps via Equation (3). Furthermore we propose an analysis, comparing classical spatial interpolation of individual GEV distributions: polynomial and spline-based regression models, inverse distance weighted and universal kriging models. This comparison is crucial to show the difference between practical commonly implemented routines and the approach proposed in this paper, using recently introduced tools, in terms of obtained return level maps (see Figure 8). The article is organized as follows. Section 2presents the considered precipitation dataset. In Section 3, return level maps are obtained through smooth spatial GEV models by using our clustering hierarchical copula-based model with associated copula-based dissimilarity measure. Obtained results on considered rainfall data are presented in Section 4. Conclusion and perspectives are discussed in Section 5. Appendix Ais devoted to illustrate, via simulation studies, the performance of the proposed dissimilarity measure and associated clustering algorithm. In Appendix Bwe provide a more comprehensive simulation study to evaluate the estimation performance for the final return level maps and compare it to alternative methods (the Gaussian latent model in Zhang et al., 2020 and the hierarchical max-stable model in Reich & Shaby, 2012). In Appendix Cwe build return level maps via classical spatial interpolation methods of individual GEV distributions. Further details on L-moments for GEV parameters are postponed in Appendix D. Finally additional information are provided in Data S1. 2CONSIDERED CENTER EASTERN CANADA DATASET We consider rainfall data measured in millimeters (mm), adjusted for snow/ice, registered in 116 stations in Center Eastern Canada for a duration of 24 h. The stations are managed by Environment and Climate Change Canada (ECCC) and verify the quality standards imposed by the World Meteorological Organization. Annual maxima precipitation for a 24-h duration are recorded from 1914 to 2017. Each station possesses precipitation measures for a minimum of 16 years in the specified time range. A specific characteristic of the considered rainfall dataset is the shortage of the recorded data, as depicted in Figure E1, illustrating the proportion of missing data and concomitant observations across years. The elevationofthestudiedareacoversawiderangegivenbetween5and672mabovesealevel.Theseannualmaximaarepublicly available at https://climate.weather.gc.ca/prods_servs/engineering_e.html. We introduce the following mathematical notation, which will be used in the following for the description of the considered models. Let xi=(x(i,1),…,x(i,T))be the annual rainfall maxima time series of the ith station, for i∈I∶= {1,…,n},observedforTyears. Then, our rainfall data can be represented as a (n×T)-matrix, X=(xi)i=1,…,n,where n=116 (number of stations) and T=104 (length of the considered whole time window). The spatial location of the considered stations is shown in Figure 1. Remark that all maps in this paper have been obtained with package leaftet in R. As suggested by Figure E1, we introduce the indicator variable Ii t={1,if x(i,t)is observed, 0,otherwise. (1) Furthermore,wedenotebyIf(resp.Iv) thenonemptysetofindices ofstationsused forthefitting(respectively,validation) procedure. Let nf=card (If)and nv=card (Iv). Obviously, If∩Iv=∅and nf+nv=n. Here, we consider nf=95 and nv=21. Note that the arbitrary choice of the nvvalidation stations can be influent in the final performance of the proposedmodels.Forthisreason,inordertotesttherobustnessoftheinvestigatedmodels,wedecidetochooserandomly 200 combinations of nf=95 and nv=21 stations to perform our study. One of the 200 considered random combinations between fitting and validation stations is displayed in Figure 2. The sensibility of the considered statistical routines with respect to the choice of combinations of nfand nvstations is investigated in Data S1. 1099095x, 2023, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/env.2795 by Universidad De Sevilla, Wiley Online Library on [19/09/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 4of30 PALACIOS-RODRIGUEZ et al. FIGURE 1 Locations of considered 116 stations in Central Eastern Canada. The latitude coordinates vary in the range [43.71,50.24] and the longitude ones in [−79.23,−54.57]. FIGURE 2 Spatial locations of one of the 200 considered combinations between fitting stations If(blue points) and validation stations Iv(red points) with nf=95 and nv=21. In the following we first consider spatial smooth GEV models to derive return level maps for every location s∈S, where Srepresents the overall surface being interpolated. The involved parameters are modeled via smooth functions in space by including some significant covariates of the models (see Section 3). Then, we compare the obtained results with classical spatial interpolation methods of individual GEV distributions (see Appendix C) based on the L-moments estimators (see Appendix D). 3RETURN LEVEL MAPS VIA SPATIALLY SMOOTH COPULA-BASED GEV MODEL In this section, we describe the components of the spatial smooth GEV model. We detail the GEV model used, the polynomial and spline models and the copula model used for the GEV parameters, which takes into account the peculiarity of the dataset presented in Section 2. 3.1 Univariate EVT via block maxima approach SinceXisdefinedbyannualmaximaofprecipitation,wefocusontheEVTblock-maximaapproach(see,e.g.,Coles,2001; Ferreira & de Haan, 2015). Let x(i,t)be the annual maxima at the ith station for year t. Then, we write x(i,t)=maxj{zj (i,t)}, 1099095x, 2023, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/env.2795 by Universidad De Sevilla, Wiley Online Library on [19/09/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License PALACIOS-RODRIGUEZ et al. 5of30 where zj (i,t)represents the precipitation at the ith station the jth day of the considered year t. EVT requires independence orshort-rangedependence(Leadbetteretal.,1983)andwedoobservethroughanadditionalanalysis,near-independence in our precipitation time-series for every considered year and station. Therefore, we can model x(i,t)by means of a GEV distribution with parameters Λ=(𝜇,𝜉,𝜎),thatis, G(x;Λ)=⎧ ⎪ ⎨ ⎪ ⎩ exp{−[1+𝜉(x−𝜇 𝜎)]−1∕𝜉},1+𝜉(x−𝜇 𝜎)>0, 0,otherwise. (2) The shape parameter 𝜉∈Rdescribing the tail behaviour of the distribution is called the extreme value index, 𝜇∈R is the location parameter and 𝜎>0 the scale parameter. Now, we introduce the return level for the GEV distribution. The return level q(p;Λ) associated with the return period 1∕p(0 <p≤1) is the (1 −p) th quantile of the GEV distribution in (2), that is, it is expected to be exceeded on average once every 1∕pyears: q(p;Λ)={𝜇−𝜎 𝜉[1−{−ln(1−p)}−𝜉],𝜉≠0, 𝜇−𝜎ln{−ln(1−p)},𝜉=0.(3) Inthefollowing,weconsiderthreeclassicalgeographicalcoordinatesascovariates:longitude,latitude,andelevation, which are obtained using a digital elevation model. Furthermore, in our analysis, we include the 75% quantile precipitation over the 34-year period 1981–2014 of the Canadian Regional Climate Model (CRCM5) driven by the ERA-Interim reanalysis (see, e.g., Dee et al., 2011). From CRM5, variables are available on a regular lattice covering the northeastern partof NorthAmerica through90,000 gridcells,each ofwhich correspondsto anareaof12 ×12km2in size.Weconsider the coordinates of the centers of these cells and we attribute to each station in Figure 1the 75%quantile precipitation value of the spatially nearest center-cell. More information on this climate reconstruction can be found, for instance, in Bresson et al. (2017). 3.2 Smooth models for GEV parameters In this section, we consider the estimation of a spatial smooth GEV distribution with the joint use of all stations. We aim to model the GEV parameters Λ(s)for s∈Sfrom the data as smooth functions in space. Let 𝜁be one of three GEV parameters(either 𝜇,𝜉or ln(𝜎))and 𝜁be the associated interpolated value. Consider the following general regression model associated to the function F 𝜁(s)=F(y(1) s,,…,y(r) s)+𝜖s,(4) where y(j) s, for j∈{1,…,r}are covariates and 𝜖sis an error term. It is assumed that 𝜖sare independent and identically distributedbyanormaldistributionwithmean0andconstantvariance.Inwhatfollows,theGEVparameter𝜁atlocation sismodeledby(4)without thestochastic(Gaussian)contributionrepresentedby𝜖s.Tocorroboratetheconstant variance hypothesis,weprovideaspatialmapof 𝜖stocheckifinourspecificdatasetthereexistssomespatialpattern.Theinterested reader is referred to Figure E2 in Data S1. 3.2.1 Polynomial regression model We consider Fin Equation (4) as linear with respect to each covariate. In the present study, we consider covariates y(j) s, for j∈{1,…,r}, as polynomials of longitude, latitude, altitude and 75% quantile precipitation with a maximum polynomial degree of 3 and we take all possible combinations between these covariates with a maximum polynomial interaction degree of 3. This provides a commonly used polynomial regression model for our predictive analysis, that is, the interpolated value at location sis written as  𝜁(s)=  𝛽0+ 𝛽1y(1) s+···+  𝛽ry(r) s,(5) 1099095x, 2023, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/env.2795 by Universidad De Sevilla, Wiley Online Library on [19/09/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 6of30 PALACIOS-RODRIGUEZ et al. with the previously described covariates y(j) sand where  𝛽0,…, 𝛽rare the classical least square estimates of regression parameters(see, e.g., Rencher & Christensen, 2012). A simplified similar framework is considered in Padoan et al. (2010) where the parameters are modeled by linear regression using latitude, longitude and elevation as covariates. 3.2.2 Spline-based regression model In order to generalize Equation (5), we can model the relation between the covariates with a smooth nonlinear function Fin (4). To avoid having to deal with the estimation of a large number of parameters, we consider here a partial linearity by the following generalized additive model: 𝜁(s)=𝛽0+𝛽1y(1) s+···+𝛽qy(q) s+F(y(q+1) s,…,y(r) s)+𝜖s,(6) where𝜖sisan errortermandFis apenalized splinewith radialbasisfunctions (see,e.g.,Marx& Eilers,1998). Therefore, the interpolated value at location sis given by  𝜁(s)=  𝛽0+ 𝛽1y(1) s+···+  𝛽qy(q) s+ F(y(q+1) s,…,y(r) s),(7) where  Fis the estimated penalized spline in Equation (6) obtained by minimizing the sum of squared errors subject to constraints on its parameters, to avoid over-fitting (see, e.g., Sec. 3 in Ruppert et al., 2003). Here we take into consideration covariate models provided by polynomial regression as in (5) or spline-based regression as in (7) with longitude, latitude, altitude and 75% quantile precipitation as covariates. Table 1gathers the considered covariate models for GEV parameters. As mentioned above, polynomial regression models in Table 1are polynomials of considered covariates with a maximum degree of 3. We take here all possible degree combinations (see details in Algorithm 5) and we select the best linear regression model with the help of AIC (see Akaike, 1974). Furthermore for the proposed spline-based regression model inTable1weconsider10,000combinations of15knotsamongnf=95fittingstationsandweselect thebestsplinemodel via the generalized cross-validation (GCV) score (see details in Algorithm 6). Table 1allows us to limit the number of possible smooth GEV parameter models by considering a total of 4 (models of 𝜇)×4(models of 𝜉)×6 (models of ln(𝜎))=96 models combinations. In Section 3.3 below, we integrate the response surfaces for modeling the GEV parameters gathered in Table 1(see also Equations 5and 7) to a flexible joint dependence framework via a hierarchical copula-based model. Although the assumption of spatial independence between the stations is very unlikely in real life, it can be found in several papers and can provide satisfying results if we fix all our interest in marginal distributions (see for instance Blanchet & Lehning, 2010). Nevertheless, the aim of Section 3.3 will be to relax in a tractable way this hypothesis to build a more realistic spatial-dependent setting. TABLE 1 Considered generalized extreme value covariate parameter models from polynomial Regression in (5) (PR model) and spline-based regression in (7) (spline model). 𝝁and 𝝃parameters ln(𝝈)parameter Considered models Considered covariates Considered covariates PR model •(long, lat, alt) •(long, lat, alt) •(long, lat, alt, 75% quantile prec) •(long, lat, alt, 75% quantile prec) •(long, lat, alt, 𝜇parameter) Spline model •(long, lat, alt) •(long, lat, alt) •(long, lat, alt, 75% quantile prec) •(long, lat, alt, 75% quantile prec) •(long, lat, alt, 𝜇parameter) Notes: The considered covariate are the three geographical coordinates (long, lat, alt), the 75% quantile precipitation (75% quantile prec) and the 𝜇parameter. 1099095x, 2023, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/env.2795 by Universidad De Sevilla, Wiley Online Library on [19/09/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License PALACIOS-RODRIGUEZ et al. 7of30 3.3 Log-likelihood for the hierarchical copula-based model In order to estimate the parameter models  𝛽jof the GEV presented in Table 1, we apply a log-likelihood approach which requiresthe jointdistributionofannualmaximumprecipitationofthe consideredfittingstationsIf.Tothisend, wefocus here on multivariate hierarchical copula models; models that are able to capture different dependencies between and withindifferentgroupsofrandomvariablesviadependencecopulafunctions.Onesuchclassofmodelsisbasedonnested Archimedeancopulas(seee.g.,Hofert &Pham,2013).A(partially)nestedArchimedeancopulaCwithtwonestinglevels and Kchild copulas(or groups), is given by C(u)=C0(C1(u1),…,CK(uK)),u=(u1,…,uK)t,(8) whereKdenotesthedimension of C0(i.e.,the number ofclusters)and eachcopulaCkis Archimedeanwith a completely monotonegenerator𝜓k,fork∈{0,…,K}(see,e.g.,Nelsen,1999).Inthefollowing,forthesakeofsimplicity,weconsider C0as the independent K-dimensional copula, that is, C0(v1,…,vK)=∏K i=1vi,withvi∈[0,1]. Definition1 (Hierarchicalcopulalog-likelihood). DenotebyGi(⋅;Λ(si)),theGEVdistributionin(2)withGEVsmooth surface parameters Λ(si)as Equations (5) and (7) associated to the ith station and gi(⋅;Λ(si))its density, for i∈If.Let Ii tas in (1). Also, let k∈{1,…,K}and Ik={i(k) 1,…,i(k) dk},wheredk=card(Ik), be the set of station indices belonging to the kth cluster, such that ∪K k=1Ik=Ifand Ik∩Ik′=∅,∀k≠k′.Denotebyc𝜃k, the Archimedean copula density for the kth cluster in the nested model in (8). Let ⊥(Λ) ∶= ∑ i∈If T ∑ t=1 s.t.Ii t=1 ln{g(x(i,t);Λ(si))}.(9) Then, we introduce the log-likelihood associated to the hierarchical copula model in (8) C(Λ) = K ∑ k=1 T ∑ t=1 s.t.Ii(k) 1 t=…=Ii(k) dk t=1 ln{c𝜃k(Gi(k) 1(x(i(k) 1,t);Λ(si(k) 1)),…,Gi(k) dk(x(i(k) dk,t);Λ(si(k) dk)))}+⊥(Λ).(10) Obviously,inthe spatialindependence settingC(Λ) in (10)reducesto⊥(Λ) in(9). Thisspatial independentsmooth GEV model ⊥(Λ) in (9) was previously proposed and discussed by Blanchet and Lehning (2010). 3.4 Adapted copula-based clustering method Partitioning Around Medoids(PAM) is a well-recognized technique to create clusters with a good partitioning using medoids for a given number of clusters K(see, e.g., Kaufman & Rousseeuw, 1987,1990, Chap. 2). The PAM algorithm is based on the search for Krepresentative objects or medoids among the observations of the dataset. After finding a set of Kmedoids, clusters are constructed by assigning each observation to the nearest medoid. Next, each selected medoid object xkand each non-medoid data point xiare swapped and the objective function is computed. The objective function usedis thesum ofan appropriatedissimilaritymeasuredik(xi,xk)computedbetweenthe timeseries ofthe ithstation and the time series of the kth medoid(Reynolds et al., 2006; Schubert & Rousseeuw, 2019). The objective is to improve the quality of the clustering by exchanging selected objects (medoids) and nonselected objects. If the objective function can be reduced by interchanging a selected object with an unselected object, then the swap is carried out. This is continued until the objective function can no longer be decreased. InthefollowingwewillrunthePAMalgorithmbyusinganadaptedversionofthecopula-baseddissimilaritymeasure recently introduced by Disegna et al. (2017) to detect clusters between spatially near and dependent stations. Using the latitude and longitude covariates, we can construct additional information on stations, constituted of an(nf×nf)data matrix S, whose generic entry sij can be interpreted as the spatial distance between the ith and jth stations and  sij =sij∕(max i,j=1,…,nfsij),(11) the normalized spatial distance between the ith and jth stations. 1099095x, 2023, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/env.2795 by Universidad De Sevilla, Wiley Online Library on [19/09/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 8of30 PALACIOS-RODRIGUEZ et al. Notice that the large proportion of missing values of our(n×T)-data matrix Xrequires some crucial adaptation of classical clustering copula methods. To this end, we now introduce the notion of the bivariate hybrid empirical copula (see e.g., Segers, 2015). Definition 2 (Hybrid empirical copula). Let iand jbe fixed, with i,j∈{1,…,nf}. Consider the 2 ×T matrix composed of (x(i,t),x(j,t))⊤ {t=1,…,T},where⊤represents the transpose operator. In each column, one or both entries may be missing. Formally, our observations consist of a sample of independent, identically distributed quadruples (Ii t,Ij t,Ii tx(i,t),Ij tx(j,t)),for t∈{1,…,T}, withIi tasin(1). Let  Λi LM the L-moments estimators of the GEVparameters relative ofthe ith station(see Appendix Dfor further details). Then, the hybrid empirical copula is defined  Cij(u,v)= H(q(u, Λi LM),q(v, Λj LM)),for (u,v)∈[0,1]2,(12) where q(⋅)is as in (3)and  HT(x,y)=∑T t=11{x(i,t)≤x,x(j,t)≤y,Ii t=Ij t=1} ∑T t=11{Ii t=Ij t=1}.(13) Thehybridempiricalcopulain(12)and(13)issimilartotheclassicalempiricalcopulaprocess,butnowtheasymptotic variances and covariances are to be multiplied by the reciprocals of the observation probabilities P[Ii t=1],P[Ij i=1]and P[Ii t=Ij t=1]. Details are given in Segers (2015). Then, the adapted empirical version of the copula-based dissimilarity measure in Disegna et al. (2017) can be defined as follows. Definition 3 (Empirical hybrid copula-based dissimilarity measure). We define the dissimilarity measure  dij =f(||𝛽(M− Cij)+ sij (1−𝛽)(M−W)||),(14) where • sij is the normalized spatial distance in (11); •Mis the Fréchet upper-bound copula, that is, M(u,v)=min(u,v); •Wis the Fréchet lower-bound copula, that is, W(u,v)=max(u+v−1,0); •𝛽∈[0,1]is the tuning parameter which reflects the prior belief of the decision maker about the desired influence of the spatial component on the clustering procedure; • Cij is the hybrid copula defined as (12) and (13); •||⋅||is a suitable norm in the copula space and fis an increasing and continuous real-valued function with f(0)=0. Note that the considered copula-based dissimilarity measure in (14) can be formalized as a suitable function of the hybrid empirical copula  Cij (expressing the dependence between the ith and jth stations) and the spatial information sij. The weight of this convex combination is expressed by the magnitude of the 𝛽 parameter. It seems that this choice f(⋅)=exp(⋅)−1 is the most convenient to highlight small differences among dissimilarity values, but clearly other functions fcould be used provided that fis increasing with f(0)=0. In Algorithm 1, we detail the first part of the proposed inference procedure, that is, the estimation of the dissimilarity measure  dij in (14) for our dataset. - First, we build the pseudo-observations of our spatiotemporal rainfall maxima, using the point-wise GEV quantiles with parameters estimated via L-moments (see steps 1 and 2 in Algorithm 1and Appendix D). - Second, we estimate the dissimilarity measure on these pseudo-observations via the hybrid copula presented in (12) and (13)(seesteps3and7inAlgorithm1). 1099095x, 2023, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/env.2795 by Universidad De Sevilla, Wiley Online Library on [19/09/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License PALACIOS-RODRIGUEZ et al. 9of30 Algorithm 1. Proposed implementation of a copula-based dissimilarity measure (Step 1) Estimate  Λi LM =(𝜇i LM,𝜎i LM, 𝜉i LM), i.e., the L-moments estimators of the GEV parameters relative of the ith station. (Step2)FromEquations(2)and(3),estimatetheinverseofthemarginalparametricestimatoroftheGEVdistribution of the daily annual maxima of the ith station, that is, q(⋅, Λi LM)and build the pseudo-observations for each station. (Step 3) Fix 𝛽∈[0,1]; (Step 4) Evaluate  sij as in (11). (Step 5) Choose the Crámer–von Mises L2norm in (14) and f(⋅)=exp(⋅)−1. (Step 6) Evaluate nc f=∑T t=11{Ii t=Ij t=1}, that is, the number of identical observations in (x(i,t),x(j,t))⊤ {t=1,…,T}. (Step 7) Fix a threshold value n. if nc f≥n, using (12) and (13) and (Step 2), evaluate  dij as in (14) on the pseudo-observations. if nc f<n,theith and jth stations are assumed to be independent. if nc f=0, instead of the dissimilarity measure dij, we only consider f( sij),thatis,fapplied on the normalised spatial distance between the ith and jth stations. The associated Rcode can be found in dissimilaritymeasure.R file in CodeR folder in Data S1. Algorithm 2. Proposed implementation for smooth copula-based GEV method (Step 1) Choose several values for K(i.e., the number of clusters) and 𝛽. (Step 2) For the input parameter 𝛽, estimate the copula-based dissimilarity measure  dij via Algorithm 1. (Step 3) Provide the clustering of stations by running the PAM algorithm for Kclusters with  dij from (Step 2). (Step 4) Estimate the Archimedean copula density c 𝜃kvia canonical maximum likelihood on the pseudo-observations for k∈{1,…,K}. (Step 5) Select the best one in terms of the AIC criterion for each cluster with k∈{1,…,K}. (Step 6) Let c 𝜃kfor k∈{1,…,K}as in (Step 5). Let define the marginal GEV smooth surfaces 𝜇(s),𝜉(s)and ln(𝜎(s)) as in Equations (5) and (7) with covariate models gathered in Table 1. Maximize the log-likelihood C(Λ) in Equation (10) with respect to the vector of parameters  𝛽in Equations (5) and (7). (Step 7) Select (K∗,𝛽∗)to minimize normalized scores in Table 2. The associated Rcode can be found in smoothDEPENDENCEmodel.R file (resp. smoothINDEPENDENCEmodel.R file for the independence setting of Equation (9)) in CodeR folder in Data S1. In Algorithm 2, we detail the second part of the proposed inference procedure, that is, the maximization of the hierarchical log-likelihood C(Λ) in Equations (9) and (10) in terms of the GEV smooth surface parameters as in (5) and (7). - For a fixed number of clusters Kand 𝛽in (14) we provide the spatial clustering of stations, using the PAM algorithm (steps 1–3 in Algorithm 2). - Then, we estimate the Archimedean copula density via the canonical maximum likelihood on the clustered pseudo-observations (see steps 4 and 5 in Algorithm 2). - Finally, we maximize the hierarchical loglikelihood (see (9) and (10)) in terms of the GEV smooth surface parameters’ via polynomial regression and spline models with the covariates presented in Table 1(see step 6 in Algorithm 2). A crucial step in the proposed construction is the optimal selection of parameters (K,𝛽)(see step 7 in Algorithm 2). To this aim we introduce in Table 2several well-known normalised scores, used in the following as selection criteria. The normalization is motivated by the spatial heterogeneity in extreme precipitation in our large region of interest. We measure via these scores the quality of the goodness-of-fit of the quantiles of the interpolated GEV parameters versus the observed ones. Let z(1) si,…,z(m) si,…,z(M) sidenote the M=30 empirical quantiles at location si, where the mth value z(m) siis associated to a probability pm=m−1∕2 M, for m=1,…,M.Also, qpm,sidenotes the (1−pm)quantile of the interpolated GEV distribution at location si, obtained by Equation (3) replacing Λby the interpolated values  Λwith p=pm. 1099095x, 2023, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/env.2795 by Universidad De Sevilla, Wiley Online Library on [19/09/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 16 of 30 PALACIOS-RODRIGUEZ et al. the spatial effect is modeled via Gaussian Markov random fields. Indeed, the 24-h duration 20-year precipitation return level map in Perreault et al. (2019) shows a very similar behavior of smooth copula-based one (left column, fourth row in Figure 8). 5CONCLUSION The statistical methodology presented in this paper, based on a hybrid hierarchical smooth GEV copula model, can be useful to infer on any location in the area of study. It takes into account the spatiotemporal structure, using the distance between dependent clusters, by means of a dissimilarity measure designed to handle missing data. Note that other papers tackle the problem of missing data, using hierarchical max-stable process construction such as in Reich and Shaby (2012) or conditioning of a latent process in the context of extremes, such as in Zhang et al. (2021). Advanced approaches for estimating the GEV parameter surfaces have recently been introduced in the literature, such as the Max-and-Smooth approach of Johannesson et al. (2021) for flood frequency data or the approach of Sass et al. (2021) based on the spatial GEV by introducing fused lasso and fused ridge penalty for parameter regularization. The hybrid hierarchical smooth GEV copula model allows the direct mapping of quantiles by dealing simultaneously with the nonconcomitant record periods between recording stations. An interesting future extension of the proposed method is the quantification of the impact of the portion of missing values and nonconcomitant record periods at different stations on the final performance of Algorithms 1and 2.This crucial point can be explored via the study of the variance of the empirical hybrid copula-based dissimilarity measure  dij in (14). Some useful preliminary results on this issue are provided by Segers (2015) and recently by Boulin et al. (2022). Then a comparison on the ability to handle missingness would be an interesting future work. ACKNOWLEDGMENTS The authors express their gratitude to two anonymous Referees and Associate Editor for their valuable comments on this article. The authors thank Jonathan Jalbert(Polytechnique Montréal, Canada) for useful discussions. Funding in support of this work was provided by the SIMONS Foundation and the Center for Mathematical Research with the Program Simons-CRM. This work has been supported by the French government, through the 3IA Côte d’Azur Investments in the Future project managed by the National Research Agency(ANR) with the reference number ANR-19-P3IA-0002. We are grateful to Mitacs Globalink Research Award for finance support. FUNDING INFORMATION Funding in partial support of this work was provided by the SIMONS Foundation and the Center for Mathematical ResearchwiththeProgramSimons-CRM.ThisworkhasbeensupportedbytheFrenchgovernment,throughthe3IACôte d’Azur Investments in the Future project managed by the National Research Agency(ANR) with the reference number ANR-19-P3IA-0002. We are grateful to Mitacs Globalink Research Award for finance support. CONFLICT OF INTEREST STATEMENT The authors have no relevant financial or nonfinancial interests to disclose. The authors have no conflicts of interest to declare that are relevant to the content of this article. All authors certify that they have no affiliations with or involvement in any organization or entity with any financial interest or nonfinancial interest in the subject matter or materials discussed in this manuscript. The authors have no financial or proprietary interests in any material discussed in this article. DATA AVAILABILITY STATEMENT The dataset is publicly available at https://climate.weather.gc.ca/prods_servs/engineering_e.html. The data that support the findings of this study are available in Climate weather Canada at https://climate.weather.gc.ca/prods_servs/ engineering_e.html. These data were derived from the following resources available in the public domain: - Climate weather Canada web, https://climate.weather.gc.ca. The code R for this work is available in CodeR folder in Data S1. 1099095x, 2023, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/env.2795 by Universidad De Sevilla, Wiley Online Library on [19/09/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License PALACIOS-RODRIGUEZ et al. 17 of 30 ORCID Fatima Palacios-Rodriguez https://orcid.org/0000-0002-8849-3019 REFERENCES Achilleos, G. A. (2011). The inverse distance weighted interpolation method and error propagation mechanism – Creating a DEM from an analogue topographical map. Journal of Spatial Science,56(2), 283–304. Akaike, H. (1974). A new look at the statistical model identification. IEEE Transactions on Automatic Control,19(6), 716–723. Beguería, S., & Vicente-Serrano, S. M. (2006). Mapping the hazard of extreme rainfall by peaks over threshold extreme value analysis and spatial regression techniques. Journal of Applied Meteorology and Climatology,45, 108–124. Beirlant,J.,Goegebeur, Y.,Segers, J., & Teugels, J.(2004). Statistics of extremes: Theoryand applications Wiley seriesin probability andstatistics. Wiley. Blanchet, J., & Davison, A. C. (2011). Spatial modeling of extreme snow depth. The Annals of Applied Statistics,5(3), 1699–1725. Blanchet, J., & Lehning, M. (2010). Mapping snow depth return levels: Smooth spatial modeling versus station interpolation. Hydrology and Earth System Sciences,14(12), 2527–2544. Boulin, A., Di Bernardino, E., Laloë, T., & Toulemonde, G. (2022). Non-parametric estimator of a multivariate madogram for missing-data and extreme value framework. Journal of Multivariate Analysis,192, 105059. Bresson, E., Laprise, R., Paquin, D., Thériault, J., & de Elía, R. (2017). Evaluating the ability of crcm5 to simulate mixed precipitation. Atmosphere-Ocean,55(2), 79–93. Burrough, P. A. (1986). Principles of geographical information systems for land resources assessment. Oxford University Press. Chilès, J., & Delfiner, P. (2009). Geostatistics: Modeling spatial uncertainty. John Wiley & Sons, Inc. Coles, S. (2001). An introduction to statistical modeling of extreme values Springer series in statistics. Springer. Cooley,D.,Nychka,D.,&Naveau,P.(2007).Bayesianspatialmodelingofextremeprecipitationreturnlevels.JournaloftheAmericanStatistical Association,102, 824–840. Cressie, N. (1993). Statistics for spatial data. John Wiley & Sons. Das, S., Zhu, D., & Yin, Y. (2020). Comparison of mapping approaches for estimating extreme precipitation of any return period at ungauged locations. Stochastic Environmental Research and Risk Assessment,34(8), 1175–1196. Davison, A. C., Padoan, S. A., & Ribatet, M. (2012). Statistical modeling of spatial extremes. Statistical Science,27(2), 161–186. deHaan,L.,&Ferreira,A.(2006).Extremevaluetheory.AnintroductionSpringerseriesinoperationsresearchandfinancialengineering.Springer. Dee,D.,Uppala,S.,Simmons,A.,Berrisford,P.,Poli,P.,Kobayashi,S.,Andrae,U.,Balmaseda,M.,Balsamo,G.,Bauer,P.,Bechtold,P.,Beljaars, A., van de Berg, L., Bidlot, J., Bormann, N., Delsol, C., Dragani, R., Fuentes, M., Geer, A., …Vitart, F. (2011). The ERA-interim reanalysis: Configuration and performance of the data assimilation system. Quarterly Journal of the Royal Meteorological Society,137, 553–597. Deheuvels,P.(1979).Lafonctiondedépendanceempiriqueetsespropriétés.untestnonparamétriqued’indépendance.Bulletinsdel’Académie Royale de Belgique,65(1), 274–292. Diggle, P., & Ribeiro, P. J. (2007). Model-based geostatistics. Springer-Verlag. Disegna, M., D’Urso, P., & Durante, F. (2017). Copula-based fuzzy clustering of spatial time series. Spatial Statistics,21, 209–225. Ferreira, A., & de Haan, L. (2015). On the block maxima method in extreme value theory: PWM estimators. The Annals of Statistics,43(1), 276–298. Gardes,L.,&Girard,S.(2010).Conditionalextremesfromheavy-taileddistributions:Anapplicationtotheestimationofextremerainfallreturn levels. Extremes,13, 177–204. Hofert, M., & Pham, D. (2013). Densities of nested Archimedean copulas. Journal of Multivariate Analysis,118, 37–52. Hosking, J. R. M., & Wallis, J. R. (1997). Regional frequency analysis: An approach based on L-moments. Cambridge University Press. Hosking, J. R. M., Wallis, J. R., & Wood, E. F. (1985). Estimation of the generalized extreme-value distribution by the method of probability-weighted moments. Technometrics,27(3), 251–261. Huser, R., Dombry, C., Ribatet, M., & Genton, M. G. (2019). Full-likelihood inference for max-stable processes. Stat,8(1), e218. Hwang, Y., Clark, M., Rajagopalan, B., & Leavesley, G. (2012). Spatial interpolation schemes of daily precipitation for hydrologic modeling. Stochastic Environmental Research and Risk Assessment,26(2), 295–320. Joe, H. (1994). Multivariate extreme-value distributions with applications to environmental data. The Canadian Journal of Statistics,22(1), 47–64. Johannesson,A.V.,Siegert,S.,Huser,R.,Bakka,H.,&Hrafnkelsson,B.(2021).ApproximateBayesianinferenceforanalysisofspatio-temporal flood frequency data. Annals of Applied Statistics,16(2), 905–935. Kaufman, L., & Rousseeuw, P. J. (1987). Clustering by means of medoids. In Y. Dodge (Ed.), Statistical data analysis based on the L1norm and related methods (pp. 405–416). Elsevier Science Ltd. Kaufman, L., & Rousseeuw, P. J. (1990). Finding groups in data: An introduction to cluster analysis. John Wiley and Sons. Khedhaouiria, D., Mailhot, A., & Favre, A.-C. (2020). Regional modeling of daily precipitation fields across the great lakes region (Canada) using the cfsr reanalysis. Stochastic Environmental Research and Risk Assessment,34(9), 1385–1405. Kohnová, S., Parajka, J., Szolgay, J., & Hlav cová, K. (2009). Mapping of Gumbel extreme value distribution parameters for estimation of design precipitation totals at ungauged sites.InBioclimatology and natural hazards (pp. 129–136). Springer. Leadbetter, M. R., Lindgren, G., & Rootzén, H. (1983). Extremes and related properties of random sequences and processes. Springer Verlag. 1099095x, 2023, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/env.2795 by Universidad De Sevilla, Wiley Online Library on [19/09/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 18 of 30 PALACIOS-RODRIGUEZ et al. Lehmann, E. A., Phatak, A., Stephenson, A., & Lau, R. (2016). Spatial modelling framework for the characterisation of rainfall extremes at different durations and under climate change. Environmetrics,27(4), 239–251. Marx, B. D., & Eilers, P. H. (1998). Direct generalized additive modeling with penalized likelihood. Computational Statistics & Data Analysis, 28(2), 193–209. Matheron,G.(1969).LeKrigeageuniversel.InFontainebleau:CahiersduCentredeMorphologieMathématique(Vol.1).ÉcoledesMinesdeParis. Nalder, I. A., & Wein, R. W. (1998). Spatial interpolation of climatic Normals: Test of a new method in the Canadian boreal forest. Agricultural and Forest Meteorology,92, 211–225. Nelsen, R. B. (1999). An introduction to copulas Lecture notes in statistics (Vol. 139). Springer-Verlag. Padoan, S. A., Ribatet, M., & Sisson, S. A. (2010). Likelihood-based inference for max-stable processes. Journal of the American Statistical Association,105, 263–277. Patton, A. J. (2006). Estimation of multivariate models for time series of possibly different lengths. Journal of Applied Econometrics,21(2), 147–173. Perreault, L., Jalbert, J., & Genest, C. (2019). Interpolation of extreme precipitation of multiple durations in eastern Canada [Conference presentation]. 11th International Conference on Extreme Value Analysis (EVA). Zagreb, Croatia. Perreault, L., Jalbert, J., & Genest, C. (2022). Interpolation of precipitation extremes on a large domain toward idf curve construction at unmonitored locations. Journal of Agricultural, Biological and Environmental Statistics,27, 461–486. Reich, B. J., & Shaby, B. A. (2012). A hierarchical max-stable spatial model for extreme precipitation. The Annals of Applied Statistics,6(4), 1430–1451. Rencher, A. C., & Christensen, W. F. (2012). Methods of multivariate analysis Wiley series in probability and statistics. Wiley. Reynolds, A., Richards, G., de la Iglesia, B., & Rayward-Smith, V. J. (2006). Clustering rules: A comparison of partitioning and hierarchical clustering algorithms. Journal of Mathematical Modelling and Algorithms,5, 475–504. Ribatet, M., Cooley, D., & Davison, A. C. (2012). Bayesian inference from composite likelihoods, with an application to spatial extremes. Statistica Sinica,22, 813–845. Rue, H., & Held, L. (2005). Gaussian Markov random fields. Theory and applications. Chapman and Hall/CRC. Ruppert, D., Wand, M. P., & Carrol, R. J. (2003). Semiparametric regression Cambridge series in statistical and probabilistic mathematics. Cambridge University Press. Rüschendorf, L. (1976). Asymptotic distributions of multivariate rank order statistics. Annals of Statistics,4(5), 912–923. Saad, C., El Adlouni, S., St-Hilaire, A., & Gachon, P. (2015). A nested multivariate copula approach to hydrometeorological simulations of springfloods:ThecaseoftheRichelieuRiver(Québec,Canada)recordflood.StochasticEnvironmentalResearchandRiskAssessment,29(1), 275–294. Sang, H., & Gelfand, A. E. (2009). Hierarchical modeling for extreme values observed over space and time. Environmental and Ecological Statistics,16, 407–426. Sang, H., & Gelfand, A. E. (2010). Continuous spatial process models for spatial extreme values. Journal of Agricultural, Biological, and Environmental Statistics,15(1), 49–65. Sass,D.,Li,B.,&Reich,B.J.(2021).Flexibleandfast spatial returnlevel estimation viaaspatiallyfused penalty. JournalofComputationaland Graphical Statistics,30(4), 1124–1142. Schlather,M. (2002). Modelsfor stationary max-stable randomfields. Extremes: Statistical Theory and ApplicationsinScience,Engineeringand Economics,5(1), 33–44. Schubert,E.,&Rousseeuw,P.J.(2019).Fasterk-Medoidsclustering:ImprovingthePAM,CLARA,andCLARANSalgorithms.InSimilaritysearch and applications (pp. 171–187). Springer International Publishing. Segers, J. (2015). Hybrid copula estimators. Journal of Statistical Planning and Inference,160, 23–34. Smith, R. L. (1990). Max-stable processes and spatial extremes [Unpublished manuscript]. Szolgay,J.,Parajka,J.,Kohnová,S.,&Hlav cová,K.(2009). Comparisonofmappingapproachesofdesignannualmaximumdailyprecipitation. Atmospheric Research,92, 289–307. Wood, S. N. (2017). Generalized additive models: An introduction with R. Chapman & Hall. Yoon, S., Kumphon, B., & Park, J.-S. (2015). Spatial modelling of extreme rainfall in Northeast Thailand. Procedia Environmental Sciences,26, 45–48 Spatial Statistics conference 2015. Zhang,L.,Shaby,B.A.,&Wadsworth,J.L.(2021).Hierarchicaltransformedscalemixturesforflexiblemodelingofspatialextremesondatasets with many locations. Journal of the American Statistical Association,117(539), 1357–1369. Zhang, Y., Tao, S., Chen, W., & Apley, D. (2020). A latent variable approach to Gaussian process modeling with qualitative and quantitative factors. Technometrics,62(3), 291–302. How to cite this article: Palacios-Rodriguez, F., Bernardino, E. D., & Mailhot, M. (2023). Smooth copula-based generalized extreme value model and spatial interpolation for extreme rainfall in Central Eastern Canada. Environmetrics,34(3), e2795. https://doi.org/10.1002/env.2795 1099095x, 2023, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/env.2795 by Universidad De Sevilla, Wiley Online Library on [19/09/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License PALACIOS-RODRIGUEZ et al. 19 of 30 APPENDIXA.SIMULATIONSTUDIES FOR  dij ANDASSOCIATEDCLUSTERINGALGORITHM In this section we illustrate the performance of Algorithms 1and 2in several simulated datasets. In this study, we focus onthecomplete datasetting,that is, withoutmissingdata. Inthiscase, Equation (13)is then equaltothe classicalempirical distribution function (see, e.g., Deheuvels, 1979; Rüschendorf, 1976). We reproduce here three possible scenarios: a situationcharacterizedbyspatialwell-distinguishedclusters(ModelM1),byspatialoverlappingclusters(ModelM2)and by a balance between Model M1 and Model M2 (Model M3). We generate K=5 clusters for models in Models M1 and M2 and K=10 in Model M3. From Definition 1, recall that Ikis the set of station indices belonging to the kth cluster with dk=card(Ik), the size of kth cluster. Furthermore here ∪K k=1Ik=80 stations and Ik∩Ik′=∅,∀k≠k′.DenotebyC𝜃kthe dk-dimensional Archimedean copula with dependence parameter 𝜃kfor the kth cluster related to model in (8). To generate the kth cluster and the associated dependent time-series of pseudo-observations of length T=40, we implement the simulation procedure gathered in the following Algorithm 3. Algorithm 3. Simulation procedure for spatial clusters and dependent time-series (Step 1) For each k∈{1,…,K}, consider dkthe size of the kth cluster. (Step 2) The two spatial coordinates mkof kth centroid are simulated uniformly in a consider spatial domain. (Step 3) Define Σkthe variance–covariance matrix relative to the kth cluster. (Step 4) The two spatial coordinates of station sk iassociated to the centroid mkare generated via bivariate gaussian vector, i.e., sk i∼2(mk,Σk), for i∈{1,…,dk}. (Step5) Generate time-series of pseudo-observations from C𝜃k(the dk-dimensional Archimedean copula associated to the kth cluster) of length T=40. A.1 Model M1: predominant spatial information In this first study, we implement Algorithm 3by considering clearly spatially separated cluster locations (see Figure A1). Furthermore, we generate the time-series of pseudo-observations by using the following copula models: −30 −20 −10 0 10 20 30 −30 −20 −10 0 10 20 30 True simulated cluster distribution FIGURE A1 Random generation of Model M1, with a predominant spatial dissimilarity between clusters. Here K=5. The cluster dimensions are d1=d2=20, d3=30, d4=d5=5. 1099095x, 2023, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/env.2795 by Universidad De Sevilla, Wiley Online Library on [19/09/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 20 of 30 PALACIOS-RODRIGUEZ et al. Cluster label Archimedean copula 𝜽kdkLocations color in Figure A1 1Frank 220Red 2Gumbel 520 Green 3 Gumbel 1.5 30 Orange 4Independent −5Blue 5 Frank 12 5 Magenta In this setting, we expect a small weight 𝛽in the convex combination in the proposed empirical copula-based dissimilarity measure in (14). This is exactly what we observe in Figure A2.LeftcolumnofFigureA2 displays the Average Silhouette Width (ASW) criterion in our PAM algorithm with dissimilarity in (14) for 𝛽=0.01 (first row) and 𝛽=0.25 (second row). For more details about definition and interpretation of ASW criterion in PAM algorithm, the reader is directed to Chap. 2 in Kaufman and Rousseeuw (1990). For each 𝛽, the associated cluster labeling of stations is represented in the right column of Figure A2. In particular, the PAM algorithm with 𝛽=0.01 perfectly identifies K=5 clusters and the correct cluster labeling of all 80 stations. As expected, when 𝛽increases the first and fourth clusters (red and blue points, respectively) are merged into an uniquecluster(seethesecondrowofFigureA2).Indeedinthiscase,theimportanceoftime-seriesdependence,increases in the convex combination. Roughly speaking, time-series exhibiting similar dependence structures, in particular in extremes, that is, Frank copula with parameter 2 and independent copula, are clustered together. A.2 Model M2: predominant dependence information In this second study, we implement Algorithm 3by considering spatially overlapping clusters (see Figure A3). Furthermore, we generate the time-series of pseudo-observations by using the same dk-dimensional Archimedean copulas C𝜃kas for the previous table of Model M1. 5101520 0.5 0.6 0.7 0.8 β=0.01 Number of clusters ASW 234567891011121314151617181920 −30 −20 −10 0 10 20 30 −30−20−100 102030 β=0.01 5 Clusters 5 101520 0.2 0.3 0.4 0.5 0.6 0.7 β=0.25 Number of clusters ASW 234567891011121314151617181920 −30 −20 −10 0 10 20 30 −30 −20 −10 0 10 20 30 β=0.25 4 Clusters FIGURE A2 Model M1 with true clusters spatial distribution as in Figure A1. Left column: ASW criterion for two values of 𝛽in (14). The true value K=5 is displayed in dashed vertical line. Right column: associated cluster labeling of the considered 80 stations. 1099095x, 2023, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/env.2795 by Universidad De Sevilla, Wiley Online Library on [19/09/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License PALACIOS-RODRIGUEZ et al. 21 of 30 −15 −10 −5 0 5 10 15 −15 −10 −5 0 5 10 15 True simulated cluster distribution FIGURE A3 Random generation of Model M2 with spatially overlapping clusters. Here K=5. The cluster dimensions are d1=d2=20, d3=30, d4=d5=5. Conversely to the previous Model M1, we expect here 𝛽≈1inEquation(14) in order to build clusters exclusively via the dependence structure criterion. Figure A4 displays the ASW criterion for 𝛽∈{0.75,0.95,0.99}(left column) and the corresponding cluster labeling of stations (right column). In particular, the PAM algorithm with 𝛽=0.99 is able to identifyK=5clustersandthecorrectclusterlabelingof90%ofstationswithrespecttothetruedistributioninFigureA3. A.3 Model M3: a balance between spatial and dependence information In this last study, we implement Algorithm 3by considering a mixture between spatial and dependence information (see Figure A5). Furthermore, we generate the time-series of pseudo-observations by using the following copula models: Cluster label Archimedean copula 𝜽kdkLocations color in Figure A5 1 Frank 2 8 Red 2Gumbel 5 8 Green 3 Gumbel 1.5 8 Orange 4Independent −8Blue 5 Frank 12 8 Magenta 6Joe 2 8 Pink 7 Clayton 10 8 Black 8 Clayton 2 8 Forestgreen 9Joe 48Purple 10 Joe 5 8 Gray For small values of 𝛽only six clusters are identified essentially by using the spatial proximity (see first row of Figure A6). By increasing 𝛽new clusters appear. In particular with 𝛽=0.25 the algorithm is able to identify two very different copula-structures: Joe copula (purple points) and Frank one (magenta points). The value 𝛽=0.5 allows us to correctly identify the number of clusters (K=10) and the correct cluster labeling of 79 stations with respect to the true distribution in Figure A5 (i.e., 98.75% of stations). If 𝛽becomes too large the algorithm produces uncorrected clusters labeling (see third row in Figure A6) or artificial small clusters (see fourth and fifth rows in Figure A6). 1099095x, 2023, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/env.2795 by Universidad De Sevilla, Wiley Online Library on [19/09/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 22 of 30 PALACIOS-RODRIGUEZ et al. 5 101520 0.25 0.30 0.35 0.40 β=0.75 Number of clusters ASW 234567891011121314151617181920 −15 −10 −5 0 5 10 15 −15 −10 −5 0 5 10 15 β=0.75 4 Clusters 5 101520 0.25 0.30 0.35 0.40 β=0.95 Number of clusters ASW 234567891011121314151617181920 −15 −10 −5 0 5 10 15 −15 −10 −5 0 5 10 15 β=0.95 4 Clusters 5 101520 0.25 0.30 0.35 0.40 β=0.99 Number of clusters ASW 234567891011121314151617181920 −15 −10 −5 0 5 10 15 −15 −10 −5 0 5 10 15 β=0.99 5 Clusters FIGURE A4 Model M2 with true clusters spatial distribution as in Figure A3. Left column: ASW criterion for several values of 𝛽in Equation(14).Thetrue valueK=5isdisplayedin dashed vertical line. Right column:associatedclusterlabelingoftheconsidered80stations. −20 −10 0 10 20 −20 −10 0 10 20 True simulated cluster distribution FIGURE A5 Random generation of Model M3 with a mixture between spatial and dependence information. Here K=10. The cluster dimensions are dk=8, for 1 ≤k≤K. 1099095x, 2023, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/env.2795 by Universidad De Sevilla, Wiley Online Library on [19/09/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License PALACIOS-RODRIGUEZ et al. 23 of 30 5101520 0.35 0.40 0.45 0.50 0.55 0.60 β=0.1 Number of clusters ASW 2 3 4 5 6 7 8 9 1011121314151617181920 −20 −10 0 10 20 −20 −10 0 10 20 β=0.1 6 Clusters 5101520 0.35 0.40 0.45 0.50 β=0.25 Number of clusters ASW 2 3 4 5 6 7 8 9 1011121314151617181920 −20 −10 0 10 20 −20 −10 0 10 20 β=0.25 7 Clusters 5101520 0.35 0.40 0.45 β=0.5 Number of clusters ASW 2 3 4 5 6 7 8 9 1011121314151617181920 −20 −10 0 10 20 −20 −10 0 10 20 β=0.5 10 Clusters 5101520 0.25 0.30 0.35 0.40 0.45 β=0.75 Number of clusters ASW 2 3 4 5 6 7 8 9 1011121314151617181920 −20 −10 0 10 20 −20 −10 0 10 20 β=0.75 10 Clusters 5101520 0.20 0.25 0.30 0.35 β=0.95 Number of clusters ASW 2 3 4 5 6 7 8 9 1011121314151617181920 −20 −10 0 10 20 −20 −10 0 10 20 β=0.95 16 Clusters FIGURE A6 Model M3 with true clusters spatial distribution as in Figure A5. Left column: ASW criterion for several values of 𝛽in Equation(14).ThetruevalueK=10isdisplayedindashedverticalline.Rightcolumn:associatedclusterlabelingoftheconsidered80stations. 1099095x, 2023, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/env.2795 by Universidad De Sevilla, Wiley Online Library on [19/09/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 24 of 30 PALACIOS-RODRIGUEZ et al. APPENDIX B. SIMULATION STUDIES FOR THE ESTIMATION OF THE FINAL RETURN LEVEL MAPS In order to evaluate the performance of the proposed spatial smooth GEV copula-model for the final return level maps, a simulation study is presented in this section. As in Appendix A, we take into account the complete data setting, that is, without missing data. Firstly, we construct the true simulated dataset where the dependent structure between stations and the spatial surfaces of GEV parameters are known. To this end, we start with the construction of the hierarchical copula model in (8)withK=7clustersof116stationslocatedasinFigureB1. Then,wegeneratetime-seriesofpseudo-observationsoflengthT=60fromthedk-dimensionalArchimedeancopulas C𝜃kas described in Table B1. Furthermore, we generate the GEV parameters surfaces by simulated realizations of Gaussian random fields and we get the simulated dataset by applying the GEV quantile in (3) to the previous time-series pseudo-observations of the previous hierarchical copula model. Using this procedure, we simulate flexible spatial surfaces for GEV parameters with nested copula spatial dependence structure. Then, we can easily calculate the true return-level maps associated to these GEV surface parameters (see Figure B4, first panel). Now,weimplementAlgorithms1and2(Steps1–3)fordifferentvaluesof𝛽andK.Particularly,wetakeK∈{1,…,7} and 𝛽∈{0.7,0.8,0.9,0.99}.InFigureB2, we present the classical ASW criterion of PAM method from Algorithm 2 (Step 3). We observe in Figure B2 that the maximum ASW value is reached for K∗=7and𝛽∗=0.7. Notice that this choice of Kand 𝛽parameters allow us to a correct cluster covering of 100% of stations with respect to the true cluster model in Figure B1. By using this clustering, we estimated by canonical maximum likelihood, the spatial copulas model and associated parameters  𝜃k. Results are gathered in Table B1 (see Algorithm 2, Step 4). FIGURE B1 Simulated clusters via the hierarchical copula model in (8) for the 116 stations. TABLE B1 Generated dk-dimensional Archimedean copulas for spatial clusters in Figure B1 and associated estimated ones (last two columns). Cluster label Color in Figure B1 True copula C𝜽kdk𝜽kEstimated copula C 𝜽k 𝜽k 1 red Joe 6 4 Joe 4.15 2blue Frank 5 5 Frank 5.58 3 yellow Gumbel 24 2 Frank 3.83 4black Frank 41 1Frank 1.16 5 magenta Joe 11 3 Joe 3.27 6orange Gumbel 22 3Gumbel 3.03 7 green Joe 7 2.5 Joe 2.47 1099095x, 2023, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/env.2795 by Universidad De Sevilla, Wiley Online Library on [19/09/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License PALACIOS-RODRIGUEZ et al. 25 of 30 FIGURE B2 AveragesilhouettewidthcriterionforthetruesimulateddatasetbyconsideringK∈{1,…,7}and𝛽∈{0.7,0.8,0.9,0.99}. TABLE B2 Normalized goodness-of-fit scores from Table 2for 21 validation stations displayed in Figure 2by using the L-moments estimators, the smooth copula-based GEV model, the Gaussian latent model in Zhang et al. (2020) and the hierarchical max-stable model in Reich and Shaby (2012). NRMSE NMAE NMPE 1. Pointwise GEV 0.07 0.04 0.12 2. Smooth copula-based model 0.28 0.14 0.64 3. Latent process 0.26 0.11 0.66 4. Hierarchical max-stable spatial model 0.37 0.22 0.71 Note: We consider here the three geographical coordinates and the 75% quantile precipitation as GEV covariate parameter models. Abbreviations: NMAE, Normalized mean absolute error; NMPE, Normalized maximum prediction error; NRMSE, Normalized root mean squared error. FIGURE B3 QQ-plots for the five validation stations plotted in Figure 2by using smooth copula-based model (green stars), the latent Gaussian model (blue stars) and the hierarchical max-stable model (red stars). The coordinates of station 1 are (47.900,−65.833), station 50 (45.283,−71.200), station 58 (45.633,−71.367), station 60 (45.500,−73.617), and station 70 (46.250,−71.217). 1099095x, 2023, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/env.2795 by Universidad De Sevilla, Wiley Online Library on [19/09/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License