Di erent nonergodic Ground-Motion Models based on spatially varying coecient models are compared for ground-motion data in Italy. The models are based di erent methodologies: Multi-source geographically weighted regression (Caramenti et al., 2020), and Bayesian hierarchical models estimated with the integrated nested Laplace approxi- mation (Rue et al., 2009). The di erent models are compared in terms of their predictive performance, their spatial coecients, and their predictions. Models that include spatial terms perform slightly better than a simple base model that includes only event and sta- tion terms, in terms of out-of sample error based on cross-validation. The Bayesian spatial models have slightly lower generalization error, which can be attributed to the fact that they can include random e ects for events and stations. The di erent methodologies give rise to di erent dependencies of the spatially varying terms on event and station locations, leading to between-model uncertainty in their predictions, which should be accommodated in a nonergodic seismic hazard assessment. ∗kuehn@ucla.edu 1 1 Introduction In recent years, there has been a push to relax the ergodic assumption (Anderson and Brune, 1999) in the development of empirical ground-motion models (GMMs). Loosely speaking, the ergodic assumption states that one can trade time for space, and that the ground-motion distribution at a particular site is the same as the ground-motion distribution averaged over many sites. This means that an ergodic GMM ignores repeatable systematic source, path, and site e ects. Due to increasing numbers of observed data, it is possible to estimate these e ects. Much has been written about nonergodic GMMs and their use in probabilistic seismic hazard analysis (PSHA). For a schematic overview of nonergodic PSHA, see e.g. Villani and Abrahamson (2015) or Walling and Abrahamson (2012). Abrahamson et al. (2019) conducted a nonergodic PSHA for three sites in California, and showed that there can be large di erences in PSHA results due to the inclusion of nonergodic e ects on ground motion. Landwehr et al. (2016) introduced the concept of spatially varying coecient models (SVCMs, Gelfand et al., 2003) for the estimation of nonergodic GMMs. In an SVCM, several coecients are functions of location, thus encoding spatial variation of ground-motion scaling. The model of Landwehr et al. (2016) is an example of a Bayesian hierarchical model, with the covariance of the random e ects modeled by a Gaussian process (GP, Rasmussen and Williams, 2006). The GP can be understood as a prior distribution over the functional form of the spatial depen- dence of the coecients, with the covariance function providing some constraints on the shape of the function. The GMM of (Landwehr et al., 2016) was based on the SVCM formulation of Bussas et al. (2015, 2017), which relied on marginalization of the random e ects. Kuehn (2021b) reformulated the model, exploiting that the number of stations and events is much smaller than the number of records. In this formulation, the size of the covariance matrix for he spatial event and site e ects depends only on the number of events and stations, while it depends on the number of records when the e ects are integrated out. This formulation is also used in recent nonergodic GMMs based on GPs (Lavrentiadis et al., 2021; Sung et al., 2021). Recently, Caramenti et al. (2020) introduced an alternative to GPs for the estimation of spatially varying GMMs, based on geographically weighted regression (GWR) (Fotheringham et al., 2002), called multi-source geographically weighted regression (MS-GWR). The method of 2 Caramenti et al. (2020) extends GWR to account for spatially varying coecients depending on both event and station coordinates. They showed an example based on peak ground acceleration (PGA) data recorded in Italy, using the ITA2018 GMM (Lanzano et al., 2019) as a reference model. Lanzano et al. (2021) investigated MS-GWR further, extending the model to di erent response spectral periods. GWR and GP based SVCMs are based on di erent underlying assumptions, which naturally raises the question how they compare. GWR constructs local models for di erent locations, based on weighting data according to their distance to the location. The weights are determined by a spatial kernel. By contrast, GP based models assume an underlying functional dependence of the spatial coecients, and place a GP prior on the functions. There have been several studies comparing Bayesian SVCMs and GWR based models (e.g. Finley, 2011; Wolf et al., 2018), which generally nd that both methods can produce similar results. None of the comparisons, though, deal with the peculiarities of empirical ground-motion data sets, with spatial processes based on di erent coordinates (event and station), and non-spatially varying random e ects. Here, I compare di erent GP based nonergodic models with each other, as well as with the MS-GWR model of (Caramenti et al., 2020). The same data set as in Caramenti et al. (2020) is used, which comprises the Italian data of the model of Lanzano et al. (2019). The GP based nonergodic models are estimated via Bayesian inference using the integrated nested Laplace approxmation (INLA, Rue et al., 2009). INLA provides an ecient way to estimate spatially varying models, and can be used to estimate nonergodic models (Kuehn, 2021b). The paper is organized as follows: I give a short overview of SVCMs and place them in the context of nonergodic GMMs. Then follows a brief introduction to the Italian data and the speci c nonergodic GMM formulation (functional form, which nonergodic terms are included). I then provide some details on the implementation of the models with INLA and MS-GWR. Finally, di erent models are compared in terms of their predictive accuracy, and their estimated spatially varying coecients. Code to run the analysis in the computer program R(R Core Team, 2021) is provided at https://github.com/nikuehn/NonErgGMM_Italy . 3 2 Spatially Varying Coecient Models In this section, I provide a brief introduction to SVCMs, and the di erences between Bayesian GP based models and GWR models. For more details, see e.g. Finley (2011); Gelfand et al. (2003) or Fotheringham et al. (2002, 2017). These models are then connected to nonergodic GMMs. A general linear model has the following form y=~ ~+ (1) whereyis the target variable, is a vector of coecients, ~is the vector of predictor variables, andis the residual. For a spatial model, where the coecients depend on location ~ x, the model becomes y(~ x) =~ 0(~ x)~(~ x) +(~ x) (2) In a Bayesian GP based SVCM, the spatial coecients are typically written as ~ 0(~ x) =~ +~ (~ x) (3) where~ (~ x) can be understood as spatial adjustment terms to the base coecients ~ . The adjustment terms are assumed to be distributed according to a GP. For computational reasons, the GPs for the di erent adjustments are often assumed to be independent (Finley, 2011), but a multivariate speci cation is possible as well (Gelfand et al., 2004). The GPs have mean zero, so the speci cation for the ith coecient becomes  iGP(0;ki(~ x;~ x0)) (4) whereki(~ x;~ x0) is the covariance function (or kernel function). Chapter 4 of Rasmussen and Williams (2006) provides an overview of covariance functions. Typically, a stationary and isotropic covariance function is employed (such as an exponential or squared exponential kernel), but nonisotropic and nonstationary models are possible (Finley, 2011; Paciorek and Schervish, 4 2006). One can understand the Bayesian SVCM as a model where the coecients are spatially correlated, and the spatial correlation structure is determined by the covariance function. The parameters of the Bayesian SVCM can be estimated via Markov Chain Monte Carlo (MCMC) sampling (Finley, 2011; Wheeler and Calder, 2007), or by optimizing the marginal likelihood (Dambon et al., 2021; Rasmussen and Williams, 2006). In GWR, the coecients for location ~ xare determined through a weighted regression, where the weights depend on the distance to the location, and are determined by a spatial kernel. Thus, the regression coecients are estimated as ~ 0(~ x) = XTW(~ x)X1XTW(~ x)~ y (5) where Xis the design matrix, ~ yis the vector of target variables, and W(~ x) is a diagonal matrix of weights. The weight matrix is location speci c, and its elements are determined by a spatial kernel, Wij=h(~ xi);~ xj. The elements of W(~ x) place greater weight on observations that are close to ~ x. Di erent kernels may be used (Fotheringham et al., 2002; Wheeler, 2014), with popular choices being the exponential or the Gaussian kernel. Apart from choosing a kernel, one must also choose a kernel bandwidth, which determines the range of in uence of the observations. This parameter is typically determined by cross-validation. In Equation (5), the bandwidth is the same for all spatial coecients. This assumption can be relaxed (Fotheringham et al., 2017). In Equation (2), all coecients ~ are spatially varying. A more exible model is one where only some of the coecients vary with ~ x y(~ x) =~ c~c(~ x) +~ 0 v~v(~ x) +(~ x) (6) where subscript cdenotes the spatially constant coecients, and subscript vdenotes the spa- tially varying coecients. Fotheringham et al. (2002) extended simple GWR to this case. In a Bayesian SVCM, this can be considered a special case of the more general formulation of Equation (2), where some of the adjustment terms are zero. 5 2.1 Nonergodic GMMs and Spatially Varying Coecient Models A general GMM has the form yes=f(~ c;~es) +Be+S2Ss+WSes (7) whereyis the ground-motion parameter of interest, f(~ c;~) is a function describing the scaling of theywith the predictor variables ~,~ cis a vector of coecients, BeandS2Ssare the event term and site term, and subscripts eandsdenote event eand station s, respectively. The event and site term are random e ects which are introduced to account for correlation of records from the same event or station. WSesis the remaining within-event/within-station residual, corresponding to in the linear model of Equation (1). It has long been recognized that the formulation of Equation (7) is overly simplistic, as it does not account for systematic, repeatable, location speci c source and path e ects. Equa- tion (7) does account for systematic site e ects via S2S, but these are speci c to each particular station, and technically do not allow to extrapolate to locations close to the station. If one accounts for systematic terms, then the GMM becomes nonergodic, which is typically written as y(~ xe;~ xs) =f(~ c;~(~ xe;~ xs)) +B 0;e+S2Ss+L2L(~ xe) +P2P(~ xe;~ xs) +WS 0;es (8) whereL2L(~ xe) andP2P(~ xe;~ xs) are systematic source and path e ects, and the dependence of the di erent variables (target, predictors, systematic terms) on event locations ~ xeand sta- tion locations ~ xsis made explicit. In this formulation, a nonergodic model consists of local adjustment terms to an ergodic base model. Often it is assumed that the event term can be decomposed into the systematic source term and the remaining event term, B=L2L+B 0 (Al-Atik et al., 2010; Lin et al., 2011). Similarly, the within-event/within-site residual is de- composed as WS =WS 0+P2P. Accounting for spatial correlation in site terms S2Scan lead to improved predictions for sites that do not have observed records (Chao et al., 2020). In an SVCM framework, the terms associated with spatially varying coecients can be 6 thought of as the nonergodic, systematic terms. The interpretation outlined above in Equa- tion (8) ts neatly with the Bayesian SVCM, as the spatially varying coecients are modeled as adjustments to a global coecient. The Bayesian SVCM is a natural extension to the simple random e ects model, which accounts for spatial e ects with a more complicated covariance structure. A GWR model is di erent, and does not lend itself to the \nonergodic adjustment to an ergodic base model" formulation. Instead, a GWR model constructs a local model for each location under consideration. One should note that neither framework is a priori better or more suited to the development of nonergodic GMMs. 3 Nonergodic Ground-Motion Models for Italy The functional form of the ITA18 model, which serves as the reference model, is as follows (Lanzano et al., 2019) log10PGA =a+b1(MWMh) 1(MwMh)+b2(MWMh) 1(Mw>Mh) + [c2+c1(MWMh)] log10q R2 JB+h2+c3q R2 JB+h2 +k log10VS30 8001(VS301500)+1500 8001(VS30>1500) +f1FSS+f2FR +B+S2S+WS (9) where 1 (MwMh)is an indicator function that evaluates to one if the condition is true, and zero otherwise. MWis the moment magnitude, RJBis the Joyner-Boore distance in km, andVS30 is the time-average shear wave velocity in the upper 30m with units of m=s.FSSandFRare indicators for strike-slip and reverse-faulting events, respectively. B,S2S, andWS are the event term, site term, and within-event/within-site residual, which are distributed according to normal distributions with mean zero and standard deviations , S2S, and SS, respectively. In their nonergodic GMM for California, Landwehr et al. (2016) made the constant coe- cient (intercept), the geometrical spreading coecient, and the anelastic attenuation coecient 7 dependent on the event location ~ xe. In addition, they had an intercept and the linear site scaling coecient that was dependent on site location ~ xs. The model of Landwehr et al. (2016) did not include an event term Bor a site term S2S. The MS-GWR model of Caramenti et al. (2020) follows the model of Landwehr et al. (2016), and includes a spatially varying geometrical spreading, anelastic attenuation, and linear site scaling coecient. It also does not include event or site terms, and does not account for spatially varying intercepts. One can rewrite the reference model of Equation (9) with spatially varying coecients as log10PGA = ~ae(~ xe) + ~as(~ xs) +b1(MWMh) 1(MwMh)+b2(MWMh) 1(Mw>Mh) + [~c2(~ xe) +c1(MWMh)] log10q R2 JB+h2+ ~c3(~ xe)q R2 JB+h2 +~k(~ xs) log10VS30 8001(VS301500)+1500 8001(VS30>1500) +f1FSS+f2FR +B 0+S2S0+WS 0 (10) where two spatial intercepts (dependent on event and station location), spatially varying geo- metrical spreading, anelastic attenuation, and linear site scaling coecients are included. If one de nes a spatially varying coecient as an adjustment to the ergodic base coe- cient,the full spatial model can be written as log10PGA =a+ae(~ xe) +as(~ xs) +b1(MWMh) 1(MwMh)+b2(MWMh) 1(Mw>Mh) + [(c2+c2(~ xe)) +c1(MWMh)] log10q R2 JB+h2+ (c3+c3(~ xe=es))q R2 JB+h2 + (k+k(~ xs)) log10VS30 8001(VS301500)+1500 8001(VS30>1500) +f1FSS+f2FR +B 0+S2S0+WS 0 (11) whereae(~ xe) is the spatially varying adjustment for the constant coecient a, dependent on 8 event location. The other adjustment coecients are denoted similarly. The subscript 0 is added to the event and site term and the within-event/within-site residual to indicate that they are the terms after the systematic, repeatable nonergodic terms are taken out. The adjustment to anelastic attenuation coecient, c3, is modeled as event location depen- dent in Landwehr et al. (2016) and Caramenti et al. (2020). This implies that the attenuation for a speci c event is isotropic, i.e. is the same in all directions. This is an improvement over an ergodic model, but is unrealistic. Hence, the anelastic attenuation coecient is made dependent on the midpoint between event and station location (denoted by ~ xes), to account for possible di- rectional di erences in attenuation. I compare the midpoint model to an event-location model, which is why the coecient is denoted as c3(~ xe=es), to specify that the adjustment depends either on event or midpoint location. Even if the anelastic attenuation is modeled dependent on the midpoint, it neglects that the anelastic attenuation is an accumulated e ect over the whole travel path. This could be modeled by a dependence on both event and station locations, but this runs into problems that there not many close points in four dimensions. Dawood and Rodriguez-Marek (2013) introduced a cell-speci c attenuation model in which the anelasic attenuation is modeled as the sum over distance fractions in small cells. Kuehn et al. (2019) cast the cell-speci c model as a Bayesian model to account for uncertainty in cells that are not well sampled. In the cell-speci c model, the adjustment to the anelastic attenuation term becomes cattn=NCX i=1cca;iRi =!R~cca (12) where Riis the fraction of the path in the ith cell. For the Bayesian SVCMs, the spatially varying adjustment coecients are modeled as a Gaussian process (or a spatial Gaussian eld, GF) with mean zero and Mat ern covariance 9 function aGP(0;!2 ak(~ x;~ x0)) (13) k(~ x;~ x0) =2(1) ()(j~ x~ x0j)K(j~ x~ x0j) (14) where!denotes the standard deviation of the GP, which is a measure of how much the coe- cient can vary across the region under study. In Equation (14),  is the gamma function, Kis the modi ed Bessel function of the second kind, is a scale parameter and is a smoothness parameter. For = 0:5, the Mat ern covariance function becomes the exponential covariance function (Rasmussen and Williams, 2006), which is extensively used in spatial correlation mod- els for ground motion (e.g. Chao et al., 2020; Esposito and Iervolino, 2011, 2012; Huang and Galasso, 2019; Jayaram and Baker, 2009; Kuehn and Abrahamson, 2020). The practical spatial range is de ned as `=p 8=, which corresponds to the distance where the correlation becomes 0.14 (Krainski et al., 2019). Here, the Mat ern covariance function with = 1 is used. This makes it possible to esti- mate the model parameters using Bayesian inference, via the INLA method (Rue et al., 2009). Together with the stochastic partial di erential equation (SPDE) approach (Bakka et al., 2018; Lindgren and Rue, 2015; Lindgren et al., 2011) INLA allows for fast and ecient inference in spatial models. Some further details on spatial models and INLA are provided in Section Im- plementation. An SVCM can be easily combined with the cell-speci c attenuation model in INLA (Kuehn, 2021b). In Equation (10), the ergodic coecients ( a,b1,b2,c1,c2,c3,k,f1,f2) are xed e ects. The adjustment terms (both the spatial adjustment terms, as well as the event terms Band site termS2S), are random e ects. The parameters that control the distribution of the random e ects are called hyperparameters. These are the standard deviations and S2Sfor event terms and site terms, and the standard deviation !and spatial range `for the spatial terms, or the standard deviation !cafor the cell-speci c adjustment coecients. All parameters ( xed e ects, random e ects, and hyperparameters) are estimated at the same time during the tting process. 10 4567 0.1 1.0 10.0 100.0 RJB (km)MW(a) Magnitude-distance distribution. 36°N38°N40°N42°N44°N46°N 6°E 8°E10°E12°E14°E16°E18°E Longitude ( °)Latitude ( °) (b) Map of data. Stations are triangles, events are cicles. Figure 1: Data used in the analysis. In the following section, I consider di erent nonergodic SVCMs, implemented in INLA. The models di er in the nonergodic adjustment terms that are included. The INLA models are compared with each other, as well as the MS-GWR model of (Caramenti et al., 2020). 4 Data The analysis is performed on the same data set as in Caramenti et al. (2020). This data set comprises the Italian subset of the data used to derive the ITA18 GMM of Lanzano et al. (2019). For details on the data set see Lanzano et al. (2019) and Caramenti et al. (2020). In total, there are 4784 records from 137 events, recorded at 923 stations. I consider stations with the same location as the same station. In the data set of (Caramenti et al., 2020), there were three stations (as identi ed by their location) with slightly di erent values of VS30. I used the geometric mean of the VS30-values of all records from each station as the nal VS30value for each station. The magnitude ranges from 4 6:87, the distance range is 0 RJB199:7km. Figure 1 shows a magnitude-distance scatterplot, as well as a map with event and station locations. For the analysis, the event and station locations are converted to UTM coordinates (UTM zone 33). The target value is horizontal PGA in units of cm=s2. The two horizontal components are combined as RotD50 (Boore, 2010). 11 5 Implementation Here, I brie y describe details on the implementation of the INLA models. All models are implemented in the Rcomputing environment (R Core Team, 2021). The implementation of the MS-GWR model follows the one provided by Caramenti et al. (2020) at https:// github.com/lucaramenti/ms-gwr/ . Code to t the INLA models can be found at https: //github.com/nikuehn/NonErgGMM_Italy . 5.1 INLA The INLA models are implemented using the R-INLA1package (Bivand et al., 2015; Lindgren and Rue, 2015; Rue et al., 2017). R-INLA is a package for approximate Bayesian inference in latent Gaussian models. For details and mathematical background on INLA, see Rue et al. (2009). For introductions to computation with R-INLA , see e.g. G omez-Rubio (2020) or Krainski et al. (2019). For a brief introduction to INLA in the context of estimation of GMMs, see Kuehn (2021a). Kuehn (2021b) showed how one can use R-INLA to eciently estimate a nonergodic GMM which is a combination of an SVCM and the cell-speci c attenuation approach. In the SPDE approach to modeling spatial data, the spatial Gaussian eld is approximated by basis functions which are evaluated on a triangular mesh (Lindgren and Rue, 2015; Lindgren et al., 2011). For a Mat ern covariance function with = 1, the eld becomes a Gaussian Markov random eld, which means that the precision matrix (the inverse of the covariance matrix) becomes sparse, which can be exploited numerically. Thus, instead of modeling the spatially varying coecients at the observation locations, they are modeled at the mesh nodes, and projected onto the observation locations. As an example, the model for the spatial adjustment term ae(~ xe) is actually ~ae(~ xmesh) =N(0;P1) ~ae(~ xe) =Ae~ae(~ xmesh) (15) 1https://www.r-inla.org/ 12 400045005000 0 500 X (km)Y (km)Figure 2: Triangular used to estimate the spatial e ects in the INLA models, generated for the Italian data, together with stations (triangles) and epicenters (circles). where P1is the covariance matrix describing the covariance of the adjustment term between the mesh nodes ~ xmesh,Pis the associated precision matrix, and Aeis a projector matrix connecting the mesh nodes to the observation (event) locations. The precision matrix Pis sparse, which can be exploited numerically (Simpson et al., 2012). Figure 2 shows the mesh used for the Italian data. For simplicity,the same mesh is used for station and event locations. The mesh has 7505 nodes, with a maximum edge length of 10km, and an outer boundary of 50km. In general, the model complexity increases with a larger number of nodes (and smaller edge lengths). Meshes with too large edge lengths will lead to biased inferences2. 2Seehttps://haakonbakkagit.github.io/btopic114.html for some informal advice on mesh generation. 13 5.1.1 Prior Distributions INLA is a Bayesian method, which means one needs to set prior distributions. For the prior distributions of the xed e ects, the following priors are used: a;b 1;b2;c1;c2;kN(0;5) c3N(0:005;0:01) f1N(0:02;0:025) f2N(0:1;0:04) The priors for the coecients associated with di erent types of faulting ( f1;f2) are based on Bommer et al. (2003). The anelastic attenuation coecient is typically small and negative, which informs the prior distribution for c3. For the other coecients, I choose a weakly infor- mative prior distribution, which means I use a prior that is not too restrictive on the parameter space. The prior distributions for the standard deviations , S2S, and SShave to be set on the log-precisions, which are internally used in INLA. A log-Gamma distribution with shape parameter 0.9 and rate parameter 0.007 is used: lnprec logGamma (0:9;0:007) lnprec S2SlogGamma (0:9;0:007) lnprec SSlogGamma (0:9;0:007) The parameters are chosen in such a way that the 5% and 95% quantiles for the prior distri- bution of the standard deviations (if the log-Gamma distribution is suitably transformed) are 0.05 and 0.45, respectively. The prior distributions for the spatial eld parameters are based on a penalized complexity 14 (PC) prior derived by Fuglstad et al. (2019) for SPDE models. The PC prior encodes a notion of parsimony, treating the spatial model as a more complex extension of a simpler base model (in this case the ergodic base model without spatial terms). The PC prior shrinks the parameters towards the simpler base model, where the spatial terms have zero marginal variance and in nite range. For details, see Franco-Villoria et al. (2019); Fuglstad et al. (2019); Simpson et al. (2017). The parameters of the PC prior for spatial models are two value uand such thatP(!>u ) = , and two values v; such thatP(`<b ) = . Thus, the PC prior is speci ed by setting our prior belief that the standard deviation of the spatial eld is larger than a certain value, and its spatial range smaller than a certain value. For the spatial range `of all spatially varying adjustment coecients, I assume the following prior P(`<30) = 0:5 meaning I assume a prior probability of 50% that the range is smaller than 30km. For the standard deviations, I assume P(!>0:3) = 0:01 for the constants ae,as, geometrical spreading adjustment c2, and linear site scaling adjust- mentk. Since the scale of the anelastic attenuation coecient is very di erent, I use a prior distribution for its standard deviation of P(!>0:01) = 0:01 For models that use a cell-speci c attenuation model, I also use a PC prior for the standard deviation!ca, which is speci ed in the same way and uses the same values. The prior distri- butions for the spatial terms are informed by results of Landwehr et al. (2016), Kuehn et al. (2019), and Kuehn (2021b), based on Californian data. 15 5.2 MS-GWR Caramenti et al. (2020) extended the simple spatial model of Equation (2) to the case where di erent coecients can depend di erent sets of coordinates (like event and station coordinates), i.e. a model like y=~ c~c+~ e(~ x)e~e(~ x)e+~ s(~ x)s~s(~ x)s (16) where subscript c,e,sdenote the sets of coecients and predictor variables that are spa- tially constant or depend on event and station coordinates. The di erent sets of coecients are estimated sequentially: the spatial coecients dependent on the station location ~ s(~ xs) are estimated rst; then, the spatial coecients dependent on the event location, ~ e(~ xe) are estimated based on partial residuals calculated from the rst step, and nally the constant coecients ~ cfrom partial residuals of the second step. Note that the order of the sequence is arbitrary, but this is the one found to be preferred by Caramenti et al. (2020), and used in Lanzano et al. (2021). For the exact derivation of the coecients, see (Caramenti et al., 2020). The event-related spatial coecients are the geometrical spreading coecient c2and the anelastic attenuation coecient c3. The spatial coecient dependent on the site location isk, the coecient controlling the VS30scaling. The spatial weights are calculated based on a Gaussian kernel. The bandwidths used for the MS-GWR model are 25km for the event-related coecients, and 75km for the station related coecient. The MS-GWR model in the speci cation of Caramenti et al. (2020) does not include sys- tematic event and site terms. Following Lanzano et al. (2021), the residuals are partitioned into between-event, station-to-station, and within-event/within-station residuals, which allows one to calculate SS, , and S2S. This model is called \MS-GWR B". The partitioning of the residuals is done using INLA. 16 Table 1: Comparison of di erent INLA models. Column Random e ects denotes which non- spatially e ect is included, column Spatial e ects which spatially varying e ect is included. A spatially varying intercept is denoted with with c, spatially vaying geometrical spreading with gs, spatially varying linear site term with vs, and spatially varying anelastic attenuation with r. Subscript e,s,esdenote dependence on event, station, and midpoint coordinate. Column Cells indicates whether cell-speci c attenuation is included. Model Random E ects Spatial E ects Cells WAIC runtime (s)† Model 1 NA NA no 3407.34 12.51 Model 2 B NA no 2398.99 34.10 Model 3 B;S 2S NA no -777.41 22.05 Model 4 B;S 2S ce, csno -1521.88 477.44 Model 5 B;S 2S ce, cs, res no -1627.29 409.91 Model 6 B;S 2S ce, cs, reno -1382.83 269.69 Model 7 B;S 2S ce, cs, gse, res, vssno 1758.34 1926.42 Model 8 B;S 2S ce, csgse, vssyes -1953.83 617.32 Model 9 B;S 2S ce, csyes -1841.30 907.45 Model 10 NA gs e, res, vssno -187.91 1199.44 MS-GWR NA gs e, re, vssno NA 11434.796 †Runtimes depend on the computational load of the computer, and can be vari- able even for di erent runs on the same computer. All rus were carrid ut on a Macbook Pro with 2.8Ghz and 16GB RAM. 6 Results 6.1 Model Comparison based on Predictive Measures First, I compare di erent (spatial and non-spatial) models t with INLA in terms of the widely applicable information criterion3(WAIC, Watanabe, 2010, 2013), which is an information criterion for Bayesian models. Similar to the Akaike information criterion (Akaike, 1973), lower values indicate a better t. The di erent models, their speci cations, and the WAIC and runtimes are listed in Table 1. The simplest (Model 1) model is a linear model of Equation (9), without any random e ects or spatially varying coecients. I then add event and site terms, as well as the spatially varying terms for the intercept (both dependent on event and station coordinates), geometrical spreading, linear site scaling, and anelastic attenuation. I also replace the spatially varying anelastic attenuation coecient with a cell-speci c attenuation term. As expected, the simplest model performs worst in terms of WAIC. Adding event terms, and in particular site terms, leads to a major improvement. Incorporating spatially varying 3also called Watanabe-Akaike information criterion 17 coecients leads to lower WAIC values, with the inclusion of the spatially varying intercepts leading to the most improvement. One can see that the model with a spatially varying linear distance term dependent on the midpoint between the event and station coordinate has lower WAIC than when the coecient depends on event coordinate. This is somewhat expected, as the anelastic attenuation should be anisotropic for a given event. The models with cell-speci c anelastic attenuation coecients perform best. This is due to the fact that the cell model treats the anelastic attenuation as a sum over the path. I also t a model that only includes spatially varying geometrical spreading, linear site scaling, and anelastic attenuation, but no intercepts or event/site terms. Such a model is closest to the MS-GWR model of Caramenti et al. (2020), which includes only these terms. In terms of WAIC, the INLA model with these terms performs better then the simplest models, but actually worse then a model without spatial terms, but site terms S2S. One should not draw conclusions about the MS-GWR model based on this comparison, though, since the WAIC values are tied to the INLA models and their parameterization. The SVCM models lead to a signi cant increase in runtime of the models. However, the runtime is still manageable, with the longest running model nishing in under 40 minutes. In fact, even the longest running INLA model compares favorable with the MS-GWR model in terms of runtime. One thing to note here is that INLA is run with a simple integration scheme and a simpler approximation to the posterior marginals4. Running Model 9 with default op- tions (which uses a composite integration scheme over the hyperparameters, and the simpli ed Laplace approximation to the posterior marginals (Martino and Riebler, 2020)) increases the runtime to 1596.59s for Model 9. The WAIC values listed in Table 1 are based on calculating the likelihood of the observations (while integrating over the posterior distribution of the parameters), adjusted by the e ective number of parameters (Vehtari et al., 2017). The likelihood of each observation is calculated conditional on all random e ects, which means the median prediction for an observation includes the event and site term, and the standard deviation is SSor 0. For ground-motion prediction, one typically wants to predict observations from new events, which means the event term is 4using control.inla = list(int.strategy = "eb", strategy = "gaussian") in the INLA call 18 0.200.250.300.350.40 Model 1Model 2Model 3Model 4Model 8Model 9Model 10MS−GWRMS−GWR B ModelRMSE(a) RMSE values for di erent models. −2000−1500−1000−5000 Model 1Model 2Model 3Model 4Model 8Model 9Model 10MS−GWRMS−GWR B Modelloglikelihood (b) Loglikelihood values for di erent models. Figure 3: Cross-validation results. For model abbreviations, see Table 1. unknown, and the aleatory standard deviation is SS=p 2 SS+ 2(or0, calculated the same way). To compare the models how they can predict new events, a 5-fold cross-validation (CV) is performed. The events are split into 5 sets. In each iteration, the observations from one set of events are left out and used as the test set, while the other sets are used as the training set. In the prediction of the test set, I include all random e ects, i.e. if there are stations with records in both the training and test, the estimated site term is used for the prediction. The root mean squared error (RMSE) is calculated over all combined test data points. Since the accurate quanti cation of the predictive uncertainty is important in PSHA, I also calculate the log-likelihood of each test observation. This indicates how well the predictive uncertainty is captured in a model, which is important for PSHA. Note that the predictive uncertainty contains both aleatory variability and epistemic uncertainty: there is aleatory vari- ability because the predictons are made for previously unseen records and evens, and there is epistemic uncertainty because the exact values of the nonergodic model terms are not known, in particular if the new data is far from observations. There is also uncertainty in the values of the xed e ects. The loglikelihood values are calculated for each observation in the test set for all 5 folds, and the total loglikelihood is computed as the sum of all observations. For the CV calculations, the coecients associated with style of faulting ( f1andf2in Equations (9) and (10)) are dropped. I only perform the CV calculations for a subset of the models of Table 1, mainly to reduce computation time. The two versions of the MS-GWR model are included. For the second model, the site terms estimated from the residuals to the 19 corresponding site terms in the test set are added to the predictions. The loglikelihood values are not calculated for the second model, since it becomes dicult to combine the uncertainty from the st model with the one from the second model. The RMSE and loglikelihood values based on CV are shown in Figure 3. The results in terms of ranking of the models are similar to the ones based on WAIC. However, one can see that the spatial models only provide a small improvement (slightly smaller RMSE, slightly larger loglikelihood) over the base model that includes event and site terms (Model 3). Based on the CV results, Model 9 performs best, which includes spatially varying intercepts and cell- speci c attenuation, but not spatially varying geometrical spreading or VS30-scaling; whereas based on WAIC, Model 8, which includes these terms, performs better. An explanation could be that an event-speci c geometrical spreading term leads to a better prediction of records within the event, but that it is not well represented by a simple spatial correlation model. Another explanation could be that there are not enough events in the data set to resolve the spatially varying geometrical spreading well enough that it proves bene cial for predicting new events. Another conclusion one can draw from Figure 3 is that while RMSE and loglikelihood lead to a similar ranking of models, there is more separation between models in the loglikelihood values. Models 4 and 8 have essentially the same RMSE, but di erent loglikelihood values, which indicates that both lead to similar median predictions, but Model 4 has a better quanti cation of the predictive uncertainty. In terms of RMSE and loglikelihood, the MS-GWR model provides an improvement over the linear model without random e ects (Model 1), but performs worse than the simple random- e ects model without spatial terms. Adding site terms to the predictions leads to improvements, with model MS-GWR B having similar RMSE values as the other models. 6.2 Standard Deviations Figure 4 compares the standard deviations of di erent models, for the event, station and record levels ( , S2S, and SS). The median of the posterior distributions are plotted, as well as the 20 0.00.10.20.30.4 Model 1 Model 3 Model 4 Model 7 Model 9Model 10 MS−GWRMS−GWR BposteriorσφSS τφS2SFigure 4: Estimated valus of the between-even standard deviation , site-to-site standard deviation S2S, and within-event/within-site standard deviation SS. For models that do not partition the residuals, the estimated total standard deviation is show. Error bars indicate the 2.5% and 97.5% quantiles of the posterior distribution. 2.5% and 97.5% quantiles. For the models that do not include event and site terms, the total standard deviation is plotted. In general, there is a reduction in the values of SS, S2S, and due to the inclusion of spatial terms, which is expected, and has been observed before (Kuehn et al., 2019; Lavrentiadis et al., 2021; Lin et al., 2011). For the MS-GWR model, the total variability is about = 0:3, which is larger than for the other spatial models. This might in part be explained that the model does not include event and site terms, which are accounted for by the spatial terms. Since the bandwidths of the spatial kernels are quite long (75km and 25km), they cannot accommodate the non-spatially varying random e ects, which would have a bandwidth (or spatial correlation length) of zero. The model that partitions the residuals of the MS-GWR model into between-event, station-to-station, and within-event residuals leads to values of that are comparable with the spatial INLA models, while SSand S2Sare larger, and are comparable with the values from Model 3, which is the non-spatial model. These results are in line with the ndings of Lanzano et al. (2021). 6.3 Spatially Varying Coecients Figures 5, 6, and 7 show maps of the spatially varying geometrical spreading, VS30-scaling, and anelastic attenuation coecients c2,k, andc3. These are shown for Model 7 (all spa- tially varying coecients), Model 8 (cell-speci c attenuation), Model 10 (no random e ects 21 or spatially varying intercepts), and the MS-GWR model. Maps for other models, as well as maps of associated uncertainties and the spatially varying intercepts, can be found at https://github.com/nikuehn/NonErgGMM_Italy . Main di erences can be seen for Model 10; since this model does not include non-spatially varying event and site terms, the model tries to accommodate these terms by very short spatial ranges for the geometrical spreading and VS30-scaling. Di erences between Model 7 and Model 8 are small for the geometrical spreading andVS30-coecient, as are di erences between the anelastic attenuation coecient for Model 7 and Model 10. Di erences between the INLA models and the MS-GWR model are more pronounced, in particular for the geometrical spreading and anelastic attenuation coecient. There are some similarities, such as low values for kin the Northwest at the border to France, and higher values in and around Tuscany. All models have positive values of the anelastic attenuation coecient in the Nothern part of Italy. This is an area where the Moho bounce is prominent (Lanzano et al., 2016), which biases the estimates of the anelastic attenuation if not properly modeled. The cell-spec c attenuation coecents (Model 8 in Figure 7) follow a similar pattern as the c3 values of Models 8 and 10, with lower values along the Tyrrhenian coast. They also show some unphysical behavior: 15 out of 603 cells have a positive cell-speci c attenuation coecient. 6.4 Prediction Here, models are compared in terms of their predictions. Predictions are calculated for an event in Middle Italy, for seismic waves traveling to the Northwest and Southeast. The geometry is shown in Figure 8. Predictions are calculated for 1km to 200km, using Model 3 (only event and site terms), Model 9 (spatially varying intercepts and cell-speci c attenuation), and the MS-GWR model. The predictions and associated epistemic uncertainties, denoted , are also shown in Figure 8. At short distances, the predictions of both SVCMs are similar, slightly higher than the ergodic base predictions from Model 3. This is due to the combined e ect of a systematic event and site term, which both models capture with their spatially varying coecients. The SVCM predictions show a di erent attenuation with distance than the base model. For the MS-GWR model, where the anelastic attenuation coecient is dependent on 22 Model 7 400042504500475050005250 −200 0 200 400 600 800 UTM X (KM)UTM Y (KM) −1.8−1.4−1.0c2(a) Model 7 Model 8 400042504500475050005250 −200 0 200 400 600 800 UTM X (KM)UTM Y (KM) (b) Model 9 Model 10 400042504500475050005250 −200 0 200 400 600 800 UTM X (KM)UTM Y (KM) (c) Model 10 (no event/site terms) MS−GWR 400042504500475050005250 −200 0 200 400 600 800 UTM X (KM)UTM Y (KM) (d) MS-GWR model. Figure 5: Values of the spatially varying geometrical spreading coecient, for di erent models. 23 Model 7 400042504500475050005250 −200 0 200 400 600 800 UTM X (KM)UTM Y (KM) −1.0−0.50.0k(a) Model 7 Model 8 400042504500475050005250 −200 0 200 400 600 800 UTM X (KM)UTM Y (KM) (b) Model 8 Model 10 400042504500475050005250 −200 0 200 400 600 800 UTM X (KM)UTM Y (KM) (c) Model 10 (no event/site terms) MS−GWR 400042504500475050005250 −200 0 200 400 600 800 UTM X (KM)UTM Y (KM) (d) MS-GWR model Figure 6: Values of the spatially varying VS30-scaling coecient, for di erent models. 24 Model 7 400042504500475050005250 −200 0 200 400 600 800 UTM X (KM)UTM Y (KM) −0.010−0.0040.002c3(a) Model 7 Model 8 400042504500475050005250 −200 0 200 400 600 UTM X (KM)UTM Y (KM) (b) Model 8, cell-speci c attenuation coe- cients. Model 10 400042504500475050005250 −200 0 200 400 600 800 UTM X (KM)UTM Y (KM) (c) Model 10 (no event/site terms) MS−GWR 400042504500475050005250 −200 0 200 400 600 800 UTM X (KM)UTM Y (KM) (d) MS-GWR model Figure 7: Values of the spatially varying anelastic attenuation coecient, for di erent models. 25 event coordinate, the predictions to the North and South are very similar; small di erences arise due to the spatially varying VS30-scaling term at the di erent stations. By contrast, the INLA model predictions diverge, due to the cell-speci c attenuation, with predictions to the South being more similar to the MS-GWR predictions. One can also see that the INLA predictions towards the South sometimes exhibit an increase in amplitudes with distance. This happens because of a positive cell-speci c attenuation coecient, which is unphysical, and should be corrected in a real application. Figure 8 also shows the standard deviation of epistemic uncertainty associated with the predictions, denoted as . Two plots are shown: in the rst one, the standard deviation associated with the systematic site terms, S2S, is included in , while in the second plot it is not. If S2Sis include in the uncertainty, then the base model (Model 3) exhibits larger values of than the spatial models. The reason is that the value of S2Sis reduced do the inclusion of the spatial terms (cf. Figure 4), while the additional uncertainty due to the spatial terms is reduced, since the predictions are made for locations close to observed data. On the other hand, if S2Sis not included, then is determined only by the uncertainty of the xed e ects, and the uncertainty due to the spatial terms. In this case, the spatial models show larger uncertainty than the base model. The values of for the INLA and MS-GWR models are similar. The value of S2Sfor the MS-GWR B model is larger than for the INLA model, which leads to a value closer to the base model when S2Sis included. In the INLA model, the uncertainty increases strongly with increasing distance. This is due to the cell-speci c attenuation model. 7 Discussion and Conclusions The main point of the present paper is to compare SVCMs for the development of nonergodic GMMs. Two methods to model the SVCMs are compared, Bayesian hierarchical models based on GPs, and local models based on GWR, all on the Italian data of (Lanzano et al., 2019). Based on the analyses carried out in this work, both methods lead to broadly similar results, but with some di erences in the details. Di erences between methods should be embraced, 26 400045005000 0 500 X (km)Y (km)(a) Mesh, data, and event/stations for which predictions are calculated. The event location is denoted by the red circle, and predictions are calculated for stations along the red and blue lines. 0123 1 10 100 RJB (km)log10 PGA Direction North South(b) Median prediction for the di erent models along the di erent directions. with φS2S 0.00.10.20.30.40.5 1 10 100 RJB (km)ψµDirection North South (c) Standard deviation of epistemic uncer- tainty associated with the predictions. in- cludes epistemic uncertainty due to spatial ef- fects, uncertainty in xed e ects, and uncer- tainty in the site terms ( S2S). without φS2S 0.00.10.20.30.40.5 1 10 100 RJB (km)ψµDirection North South(d) Standard deviation of epistemic un- certainty associated with the predictions.  does not include S2S. Figure 8: Median predictions and associated epistemic uncertainty for di erent models. The solid lines are for INLA Model 9, the dashed lines are for the MS-GWR model. The black line corresponds to Model 3 (the non-spatial base model). 27 as none of the models are a \true" representation of reality. Epistemic uncertainty, both on model parameters as well as on di erent models, is typically included in a standard PSHA. For nonergodic GMMs, the importance of accounting for epistemic uncertainty in the nonergodic terms has long been recognized (e.g. Walling and Abrahamson, 2012). Nonergodic GMMs typically provide a method to calculate within-model uncertainty due to the spatial e ects (e.g. Caramenti et al., 2020; Landwehr et al., 2016; Lavrentiadis et al., 2021), which can be incorporated into a nonergodic PSHA (Abrahamson et al., 2019). However, between-model uncertainty is generally not Both GWR and GP based SVCM methods have their strengths and weaknesses. Caramenti et al. (2020) and Lanzano et al. (2021) argue that GWR as a method is easier to understand. This is true, as in its simplest form it is a simple weighted linear regression, and the extension to MS-GWR, which takes into account multiple sets of coordinates, is relatively straightforward. A Bayesian SVCM, on the other hand, places a GP prior on the functional dependence of the spatial adjustment on the locations. However, one can understand the GP as a framework to extend the simple random e ects structure of a typical GMM to a more general structure that takes more complicated correlations between records into account. Random e ects model are common and familiar in GMM development, so this perspective should make the model easier to understand. Results based on CV have shown that inclusion of site terms S2Sis important, with Model 3 (a non-spatial model with site terms) having similar performance as the spatial models. The GP based models can include random e ects such as event and site terms directly, while they need to be estimated from partitioned residuals for MS-GWR. Figure 4 shows that the values of the resulting standard deviations SS, , and S2Sare di erent. This could be attributed to just a case of between-model uncertainty. However, the Bayesian models show a reduction in all components of the standard deviation ( SS, S2S, and ), while for the MS- GWR model the main reduction is in the value of . Given that the spatial models are modeling nonergodic source, path, and site e ects, one would expect a decrease in all components (Lin et al., 2011). Even if one accepts that nonergodic adjustment terms (and SVCMs) are just a simple extension to a traditional random e ects model, GP based models are not easy to t. Often, 28 these models are t with MCMC sampling (Finley, 2011). Since GPs scale poorly to large data sets, approximations need to be made with increasing number of observations, such as a nearest neighbourhood GP (Datta et al., 2016; Finley et al., 2018, NNGP, ). Here, I have used INLA, which allows to t the spatial models eciently. Neither INLA nor MCMC are as easy to apply as a simple (weighted) linear regression, and require some experience. However, this may be o set by the rich inferential framework supplied by the Bayesian methods, in particular the possibility to model complicated random e ects structures. This includes the cell-speci c attenuation term, which can be formulated as a random e ect, is a more physical representation of the path e ect than a simple spatially varying coecient. In addition, one can set prior distributions on all parameters, which can help constrain the model. Informative prior distributions on model coecients have been proven to be successful to improve model performance (Kowsari et al., 2020). In terms of predictive performance, the spatial models estimated with INLA outperform the MS-GWR models. Once site terms are accounted for in model MS-GWR B, the di erence in terms of RMSE based on CV is small. However, the di erent models exhibit di erent spatial patterns of the spatially varying coecients (cf. Figures 5, 6, 7), which can lead to quite di erent predictions (Figure 8), in particular at large distances. Given the similarity in RMSE, one should probably consider the notion of using di erent models based on di erent methodologies in a PSHA. Both GWR and Bayesian SVCMs employed in this work are relatively simple. The MS- GWR model uses a single kernel, and the same bandwidth for the spatially varyig coecients that depend on the same location (geometrical spreading and anelastic attenuation). Murakami et al. (2019) note \the importance of scale in spatially Varying Coecient Modeling" and note that the inclusion of di erent bandwidths for di erent coecients is generally bene cial. This can be incorporated into a GWR model based on the methodology of Fotheringham et al. (2017). In the INLA models, di erent values of `(the practical spatial range) are estimated for the di erent spatial e ects, which adds exibility to the model. Other extensions to increase exibility are the use of adaptive kernels (where bandwidth is not measured in term of length, but n the number of nearest neighbours (Wolf et al., 2018)) for GWR, or the use of non- 29 stationary and non-isotropic covariance function in Bayesian SVCMs. Apart from a comparison of di erent SVCMs, this paper also highlights the usefulness of INLA to estimate nonergodic GMMs. A range of di erent models were t with INLA, with di erent numbers of spatial random e ects. INLA is quite fast (cf. Table 1), which is important for scaling the method to larger data sets. The Italian data set used here is relatively small (4784 records) compared with other regions, which means that model tting can be done relatively easy. For comparison, the Californian data of the NGA-W2 data set Ancheta et al. (2014); Bozorgnia et al. (2014) comprises more than 12000 records. A nonergodic model for California similar to the ones estimated here takes about 1-2 hours to t with INLA, compared to 20 hours using MCMC sampling (Kuehn, 2021b). Data sets are going to increase in the next years, both due newly recorded data, as well as lowering magnitude thresholds. This is an opportunity to develop better nonergodic models (more data mans that the nonergodic terms can be better constrained), but this requires estimation methods that can handle large spatial data sets. INLA can t a nonergodic GMM for an induced data set in Oklahoma, which contains more than 90000 records from about 3600 earthquakes (Rennolet et al., 2018), in about 3 hours (M. Walling, personal communication). Thus, I believe INLA is a promising tool for the development of nonergodic GMMs for large data sets. Code to t the INLA models presented herein is available at https://github.com/nikuehn/NonErgGMM_Italy . Acknowledgements CV/MSLL and loglikelihood calculation follows ideas from the INLA Google group ( https:// groups.google.com/g/r-inla-discussion-group/c/K49zsVxOYiw/m/GwEQGk1dAwAJ ). Plots are made using the Rpacakge ggplot2 (Wickham, 2016). The mesh plot was made via the package inlabru (Bachl et al., 2019). Other packages used are sp(Bivand et al., 2013), sf (Pebesma, 2018), rgdal (Bivand et al., 2021). The MS-GWR code uses internally the package GWmodel (Gollini et al., 2015). Spatial plots (Figures 5, 7, and 6) follow the style of Caramenti et al. (2020), from https://github.com/lucaramenti/ms-gwr/ . I would like to thank Melanie Walling, Greg Lavrentiadis, Yousef Bozorgnia, Xiaofeng Ming, 30 and Elnaz Seylabi, and the nonergodic working group at UC Berkeley/UCLA for fruitful dis- cussions about INLA and nonergodic GMMs. Partial support of Paci c Gas & Electric Company and California Department of Trans- portation are gratefully acknowledged. Any opinions, ndings, and conclusions or recommen- dations expressed in this material are those of the author and do not necessarily re ect those of the sponsors.