A Bayesian method to include spatial correlation structure of residuals in empirical ground-motion models is presented, based on the integrated nested Laplace approxima- tion. The method is evaluated on a simulated data set as well as Italian strong-motion data. Results based on the synthetic data set indicate that the parameters of the model, as well as associated uncertainties, are well captured. The presented method facilitates the use of more complex models, such as accounting for nonergodic source, path, and site effects, as well as non-stationary correlation functions. Excluding spatial correlations in the empirical GMM leads to an increased value of the between-event variability; however, the total variability tends to be the same. GMMs that include spatial correlation lead to better predictive performance for within-event predictions. On the other hand, when predicting for new events, performance decreases slightly compared to models that do not account for spatial correlations. 1 Introduction Empirical ground-motion models (GMMs) are typically estimated via regression from observed strong-motion data sets. Modern GMMs account for systematic event and site terms by mod- eling these as random effects, sometimes also including regional random effects to account for regional variation in ground-motion scaling ( Stafford , 2014). While it is generally acknowledged that ground motions observed at different stations are spatially correlated, a spatial correlation structure is usually not included during the regression stage. Spatial correlation models are however estimated from residuals of a GMM (e.g. Esposito and Iervolino , 2011; Heresi et al. , 2022; Jayaram and Baker , 2009). Such models are useful for e.g. regional risk analyses ( Jayaram and Baker , 2010a; Manzour et al. , 2016; Sokolov and Wenzel , 2011) and shakemaps ( Sgobba et al. , 2023; Worden et al. , 2010). In general, a model that accounts for spatial correlation should be preferable, since it adheres closer to reality. ∗kuehn@ucla.edu , ORCID 0000-0002-3512-5300 1 Jayaram and Baker (2010b) proposed an iterative multi-stage method to include spatial correlations in mixed-effects regression models. Ming et al. (2019) found that the multi-stage method can lead to inconsistent results and is sensitive to initial, and proposed a scoring algorithm to overcome these shortcomings. Here, I show how one can incorporate a spatial correlation structure of residuals into a Bayesian GMM estimated with the integrated nested Laplace approximation (INLA, Rue et al. , 2009), which is a method for approximate Bayesian inference for latent Gaussian models. The INLA methodology has been applied in multiple fields, e.g. epidemiology ( D’Angelo et al. , 2021), air pollution mapping ( Lu et al. , 2023), ecology (Fichera et al. , 2023), modeling of seismicity ( Bayliss et al. , 2020; D’Angelo et al. , 2020; Naylor et al. , 2022) and seismic inversion ( Zhang et al. , 2016). For applications in ground-motion modeling, see e.g. Kuehn (2023); Macedo and Liu (2022); Sung and Abrahamson (2022); Walling et al. (2021). In particular, INLA easily allows to incorporate multiple random effects as well as spatial structures such as spatially varying coefficient models, which have been used to estimate nonergodic GMMs ( Landwehr et al. , 2016; Lavrentiadis et al. , 2022). Since INLA is a Bayesian method, it allows to incorporate prior information via the prior distributions for the parameters of the model, which can help to stabilize inference. In addition, priors can help to constrain parameters to lead to reasonable extrapolation. The paper is organized as follows: First, I give a brief overview of GMMs, as well as how spatial correlation can be incorporated. Then, the data and underlying model used in this study are described, both based on the model of Lanzano et al. (2019). A simulation study follows, showing that the implementation can recover the parameters in a controlled environment. Finally, the INLA spatial correlation GMM is applied to the real data. 2 Spatial Correlation Ground-Motion Model In general, a GMM has the following form Yes=f(⃗ c;⃗ xes) +ϵes (1) where Yesis the ground-motion parameter of interest (e.g. logarithmic peak ground acceleration (PGA)), eandsare indices for event and station, respectively, f(⃗ cT;⃗ xT es) is a function with coefficients ⃗ crelating the predictor variables ⃗ xtoY, and ϵis the residual. To account for cor- relation between records from the same event or station, the residual ϵis typically decomposed into event terms δB, site terms δSand within-event/within-site residuals δWS (Al-Atik et al. , 2010) ϵes=δBe+δSs+δWS es (2) Each of these terms is assumed to be normally distributed with mean zero and standard devi- ation τ,ϕS2S, and ϕSS: δB∼N(0, τ) δS∼N(0, ϕS2S) δWS∼N(0, ϕSS)(3) Often, the residual is only decomposed into event terms and within-event residuals, without site terms. In this case, the within-event residual is called δW, and the associated within-event standard deviation is ϕ=p ϕ2 S2S+ϕ2 SS. 2 The systematic event and site terms are often called random effects to distinguish them from the coefficients, which are also called fixed effects This distinction does not really make sense in a Bayesian context, since all parameters to be estimated are treated as random variables in a Bayesian model. In the following, the systematic adjustment terms are called random or latent effects. Parameters that control the distribution of another parameter, like standard deviation τfor event terms δB, are called hyperparameters. The hyperparameters also include the parameters of the spatial correlation structure, introduced later. If the latent effects are integrated out, then the model for all records can be written as a multivariate normal distribution (MVN) ⃗Y= MVN( ⃗f(⃗ x;x),Σ) Σij=δe(i),e(j)τ2+δs(i),s(j)ϕ2 S2S+δijϕSS(4) where ⃗Yis the vector of target variables, ⃗f(⃗ x;x) is the functional form, evaluated over the matrix of predictors (one row per record), and Σis the covariance matrix. The entries of the covariance matrix are non-zero for records from the same station or event. In Equation (4), δij is the Kronecker delta, which is one for i=jand zero otherwise. The index eiconnects record ito event e; similarly for stations. Including spatial correlations of of observations adds a term in the covariance matrix for records from an event that depends on the distance between the recording stations ( Ming et al. , 2019) Σij=δe(i),e(j)τ2+δs(i),s(j)ϕ2 S2S+δijϕ0+δe(i),e(j)ϕ2 ck(⃗ts(i),⃗ts(j)) (5) where ⃗ts(i)is the coordinate of the station corresponding to the ith record, and k(⃗ts(i),⃗ts(j)) is a correlation function. The inclusion of a new term means that the variance of the within- event/within-site residuals is reduced, which is now called ϕ2 0; consequently, the remaining residual is called δWS 0. If no site terms are considered, then the covariance matrix becomes block-diagonal ( Abrahamson and Youngs , 1992) (assuming that the data set is ordered by events), and the remaining residual is δW 0if spatial correlation is included, and δWotherwise. The correlation function determines the spatial correlation between two records, and typ- ically depends on the distance between the stations. If the correlation function only depends on the distance ⃗ti−⃗tj , then it is called stationary. One can rewrite Equation (1) to include all latent effects as Yes=f(⃗ c;⃗ xes) +δBe+δSs+δCes+δW 0,es (6) where δCesis a random affect associated with a zero-mean latent spatial Gaussian process (GP, also called a Gaussian field GF Rasmussen and Williams , 2006) that models the spatial correlations between records δC∼GP(0, ϕ2 ck(⃗t,⃗t′)) (7) Anderson and Uchiyama (2011) and Kuehn and Abrahamson (2020) have proposed to use spatial correlations of within-event/within-site residuals to account for path effects. Dawood and Rodriguez-Marek (2013) proposed a cell-specific attenuation model to model path effects, where the anelastic attenuation is modeled as the sum over discrete path segments through small cells, each of which is associated with its own anelastic attenuation coefficient. Such a model can be implemented as a random effects model, with an average attenuation coefficient and 3 cell-specific adjustments ( Kuehn , 2023; Kuehn et al. , 2019), which are assumed to be normally distributed with mean zero and standard deviation σcell. In the following, regression models including different random effects structures (with and without spatial correlations, with and without cell-specific attenuation) are implemented and compared. In general, a spatial correlation model should correspond to the partition of residuals used in the estimation of a GMM, which typically include site terms for modern GMMs. Most published spatial correlation models are estimated for within-event residuals δW, i.e. neglect site terms. Stafford et al. (2019) provides an example of a spatial correlation model that accounts for the full partition of residuals used in their GMM. Jayaram and Baker (2009) find that the spatial correlation is region-to-region dependent based on local geology, depending on the spatial correlation of VS30(time-averaged shear wave velocity in the upper 30m). Accounting for site terms should take this effect out of the residuals. 2.1 Implementation in INLA Spatial models of the form like Equation (6) can be easily estimated using INLA ( Rue et al. , 2009, 2017). For mathematical details of INLA, see ( Rue et al. , 2009) or Martino and Riebler (2020). Spatial modeling in INLA is based on the stochastic partial differential equation (SPDE) approach ( Bakka et al. , 2018; Lindgren and Rue , 2015; Lindgren et al. , 2011). See Krainski et al. (2019) for an introduction to spatial modeling with SPDEs and INLA. Kuehn (2023) and Lavrentiadis et al. (2022) provide INLA code for GMMs in a nonergodic context, with spatially correlated event and site terms as well as cell-specific attenuation. The spatial Gaussian field δCis approximated by basis functions evaluated on a triangular mesh ( Lindgren et al. , 2011), which are then projected onto the observed locations ⃗δC(⃗tmesh) =N(0,P−1) ⃗δC(⃗t) =A⃗δC(⃗tmesh)(8) where ⃗tmeshare the locations of the mesh nodes, Ais a projector matrix connecting the mesh nodes to the observations, and Pis the precision matrix (the inverse of the covariance describing the correlation between the mesh nodes). The correlation function k(⃗t,⃗t′) is a Mat´ ern correlation function k(⃗t,⃗t′) =2(1−ν) Γ(ν)(κ|⃗t−⃗t′|)νKν(κ|⃗t−⃗t′|) (9) where Γ is the gamma function, Kνis the modified Bessel function of the second kind, κis a scale parameter and νis a smoothness parameter. The default value in INLA for the smoothness parameter is ν= 1 for a two-dimensional field. Most ground-motion spatial correlation models use the exponential correlation function k(⃗t,⃗t′) = exph −|⃗t−⃗t′| ℓi , which is a special case of the Mat´ ern function for ν= 0.5. A helpful quantity is the practical spatial range ℓ=√ 8ν/κ, which is the distance where the value of the correlation function is about 0 .139. Note that typically studies of spatial correlations of ground motion describe the spatial correlation in terms of the effective range, which is the distance where the correlation reaches a value of 0 .05. In the following, I generally use the practical spatial range, as this is how the length-scale of the spatial field is parameterized in INLA. 4 All models are implemented in the computer environment R(R Core Team , 2021) with package R-INLA1. See G´ omez-Rubio (2020, chapter 3) for more information on mixed-effects models in INLA, and Kuehn (2021b) for an introduction to simple GMM modeling with INLA. The basic functional form is defined as a formula of linear predictors, and the event and site terms are defined as iidrandom effects, with event and station ID as the grouping factor. The spatial term δCis defined as an spde random effect. To model spatial correlations of within-event residuals δWor within-event/within-site residuals δWS , a separate spatial field ⃗δCfor each event is needed. This is achieved using either replicate orgroup (with an iid group model) in the specification for the spatial random effects. Since an individual latent spatial field is estimated over the mesh for each event, ideally the extent of mesh small is kept small. Thus, instead of working with actual site coordinates, the site coordinates relative to the event epicenter are used, i.e. the origin ( ⃗t= [0,0]T) for each event coordinate system is the epicenter. This does not affect the distances between sites, which is what is needed to calculate the correlation using the Mat´ ern correlation function. 3 Ground-Motion Model and Data The INLA regression model with spatial correlations of residuals is illustrated on the data set used by Caramenti et al. (2022), which consists of the Italian subset of the data used in the ITA18 GMM ( Lanzano et al. , 2019, 2022). The functional form of the ITA18 model is used to model the dependence of the target variable on the predictor variables. For details on model and data, see Caramenti et al. (2022) and Lanzano et al. (2019). The ITA18 GMM depends on moment magnitude MW, Joyner-Boore distance RJB, faulting style, and VS30. The functional form of ITA18 relating the target variable to the predictors is f(⃗ c;⃗ x) =a+b1(MW−Mh) 1(Mw≤Mh)+b2(MW−Mh) 1(Mw>Mh) + [c2+c1(MW−Mref)] log10q R2 JB+h2+c3q R2 JB+h2 +k log10VS30 8001(VS30≤1500)+1500 8001(VS30>1500) +f1FSS+f2FRV(10) with⃗ c= [a, b 1, b2, c1, c3, k, f 1, f2]Tand⃗ x= [Mw, RJB, VS30].FSSandFRVare indicator variables for strike-slip and reverse faulting, respectively. The parameters Mh,Mref, and hare fixed, making the model linear in the predictors. The ITA18 model is developed is developed for PGA and spectral acceleration for 26 periods from T= 0.01stoT= 10s. In this work, the regression models including spatial correlations are carried out for PGA. Figure 1 shows a magnitude-distance scatterplot, a map of the events and stations, as well as the mesh used for the SPDE approximation. In total, the data set consists of 4784 records from 137 events and 923 stations. The mesh has 4929 nodes and 9728 triangles. For the prior distributions of the coefficients, the INLA default priors are used, which is a normal distribution with mean zero and precision zero for the intercept a, and precision 0.001 for the other coefficients. In INLA, variabilities for the likelihood and iid random effects are internally represented as logarithmic precisions (which is one over the variance), so priors 1https://www.r-inla.org/ 5 4567 0.1 1.0 10.0 100.0 RJB (km)MW 36°N38°N40°N42°N44°N46°N 6°E 8°E10°E12°E14°E16°E18°E LongitudeLatitude −200−1000100200 −200 −100 0100 200 X (km)Y (km)Figure 1: Left: Magnitude distance distribution of data set. Middle: Map of events (circles) and stations (triangles). Right: Mesh used for the spatial field, based on centered station coordinates. Points are stations, color-coded by event. need to be set on this scale. The prior for the log-precision of the likelihood term (i.e. the within-event/within-site residuals) is a the default prior, a log-Gamma distribution with shape parameter 1 and rate parameter 10−5. The prior distribution for the logarithmic precisions of the event and site terms is a log-Gamma distribution with shape 1 .74 and rate 0 .0153, which is chosen such that the prior 1% and 99% quantiles for the corresponding standard deviations τandϕS2Sare 0 .05 and 0 .4. For the prior distribution of the parameters of the spatial field a penalized complexity (PC) prior ( Fuglstad et al. , 2019; Simpson et al. , 2017) is used with the following parameterization p(ℓ <50) = 0 .5 p(ϕc>0.4) = 0 .01(11) For the cell-specific attenuation model, the path length within each cell is divided by 100, to avoid very small cell-specific adjustment coefficients. The prior distribution for the standard deviation of the cell-specific attenuation coefficients is a PC prior with p(σcell>1) = 0 .01 (12) 4 Results using a Synthetic Data Set In this section, the model is evaluated on a simulated data set, for which coefficients, latent parameters, and hyperparameters are known, and the data generating process is exactly the same as in the regression the model. In that way, one can assess whether the model can recover the parameter values for the given data set, as well as what happens under model- misspecification. I take the Italian data of the ITA18 model, fix the coefficients and hyperparameters (stan- dard deviations and spatial range) and randomly sample observations. Hence, I calculate median predictions for each record, and sample event terms, site terms, spatially correlated observations, and remaining residuals, and combine them all into one data set. The coefficients are fixed to the values of ITA18, and Table 1 shows the values used for the hyperparameters. Then, regressions are performed according to the models presented previously. Different mod- els are estimated: (1) a standard GMM (ignoring spatial correlation) which only accounts for event and site terms; (2) the full model. The spatial correlation parameters are also estimated from the within-event/within-site residuals. For comparison, models that do not account for site terms are also investigated. 6 Table 1: Values of hyperparameters for simulation. τ ϕ S2SϕSS ϕcℓ 0.17 0.2 0.1 0.22 40 Table 2: Fraction of estimated coefficients within 5% to 95% quantiles of posterior distribution, over 80 different simulations, for the base model without spatial correlations, and the full model including spatial correlations. Coefficient Base Model Full Model a 0.65 0.83 b1 0.53 0.78 b2 0.85 0.92 c2 0.62 0.87 c1 0.38 0.85 c3 0.73 0.93 f1 0.88 0.82 f2 0.92 0.88 k 0.87 0.90 Table 2 shows results with respect to the coverage of the posterior distribution of the coefficients. 80 different data sets are simulated, and regressions are carried out using the base model (without spatial correlations) and the full model including spatial correlation. I then check whether the true coefficient is inside the 90% credible interval, i.e. is within the range between the 5% and 95% quantile of the (marginal) posterior distribution. In general, the model can recover the coefficients well, and the coverage for the full model is reasonably well matched with expectations, while it is generally lower for the base model. In the subsequent paragraphs, the estimated hyperparameters for one specific simulated data set are investigated. Figure 2 shows posterior distributions of the estimated standard deviations ϕ0,ϕS2S, and τfor the full model, and ϕSS=p ϕ2 0+ϕ2 cfor the standard model without spatial correlation. The true values are well within the body of the posterior distributions, indicating that the model can recuperate the parameters well (in the controlled simulation environment). When not accounting for spatial correlation, the overall variability is well captured, but the between- event standard deviation is overestimated, while ϕSSis underestimated. This effect is also observed by ( Jayaram and Baker , 2010b) and Ming et al. (2019). The underestimation of τ when not accounting for spatial correlation is due to the fact that the covariance matrix for the records from one event is Σij=τ2+ϕ2 ck(⃗ti,⃗tj), so the spatial correlation structure increases the entries. This effect increases for longer length scales, since these generally increase the value of the correlation, which is shown in Figure 3. Here, the estimated values of the standard deviations are shown (averaged over 10 simulations), for regressions that do not account for spatial correlation and for simulations with different spatial ranges. The ratio of the estimated to true variance τ2increases to almost a factor of two for longer spatial ranges, while ϕSS decreases. The single-station standard deviation σSS=p τ2+ϕ2 SSis well estimated (the increase in τand decrease in ϕSScancel out), as is the station-to-station variability ϕS2S. Figure 4 shows the posterior distribution of the spatial range ℓand associated standard deviation ϕc, for the full model, and models estimated from within-event residuals δWand within-event/within-site residuals δWS . For the latter two models, first a standard regression 7 50100150 0.10 0.15 0.20 standard deviationposterior densityφ0 φS2SσSSτ full no corrFigure 2: Posterior distributions for standard deviations τ,ϕ0,ϕS2S, and ϕSSfrom regression on simulated data, with and without accounting for spatial correlation. Vertical lines are true values used in simulation. 0.160.200.240.28 40 60 80 100 120 140 spatial range (km)standard deviationφS2S φSSσSSτ 0.81.01.21.41.61.8 40 60 80 100 120 140 spatial range (km)variance ratioφS2S φSSσSSτ Figure 3: Estimated value of standard deviations when neglecting spatial correlation, for dif- ferent values of he spatial range used in the simulation. Horizontal lines are true values. Right: ratio of variance between estimated and true values. Shown are results averaged over 10 simu- lations. 8 0.050.100.150.200.25 25 30 35 40 45 50 spatial range (km)posterior densityδW δWS full model 306090 0.17 0.19 0.21 0.23 0.25 φcposterior densityδW δWS full modelFigure 4: Posterior distributions for practical spatial range and standard deviation ϕc, estimated using the full model and based on δWandδWS from a standard regression. Vertical lines show the true values used in the simulation. 0.00.51.01.52.0 10 100 1000 10000 spatial range (km)densityinla variogram 110100 10 100 number of recordsspatial range (km) Figure 5: Kernel densities of even-specific range estimates for simulated data, estimated on δWS using INLA and variograms. Right: Range dependent on number of records per event. True value of range is 40km and is shown as black vertical/horizontal line. For INLA, the mean of the posterior distribution is used. without spatial correlation is performed, and then the spatial correlation parameters are esti- mated from the residuals. Using δWS , which are the residuals which are spatially correlated in the simulation, leads to a reasonable estimation of the spatial range, but underestimates the variability. This is due to the fact that the spatial correlation structure is estimated on point estimates of the residuals, which ignores uncertainty in the event terms δBand site terms δS, which can be quite large in particular for sites that do not have many recordings. Estimation of the spatial range using within-event residuals δWleads to a smaller value than the one used in the simulation. However, in this case the spatial range is measuring a different signal. The site terms, which are not spatially correlated in the simulation and are not considered in the estimation, add some uncorrelated noise to the residuals, which leads to some dilution of the signal and hence an underestimation of the spatial correlation. Hence, for δWthe spatial range is more like an intermediate of the spatial correlation of δWS andδS. For Figure 5, the spatial range is estimated separately for each event with at least 10 record- ings, and the within-event/within-station residuals δWS . The ranges are estimated both using INLA, and using the traditional geostatistical approach of a variogram ( Oliver and Webster , 2014). The histograms of the estimated ranges are shown in Figure 5. On average, the true value is captured quite well, though there is some variability in the estimates (recall that the true value is 40km for all events). The INLA estimates show a little bit lower variability, which is probably due to a regularization effect of the prior distribution. There is a dependence in the variability on the number of records per event, in, but even the most recorded events show some variability around the true value. 9 3.33.43.53.63.7 M0 M1 M2 M3a 0.10.20.30.4 M0 M1 M2 M3b1 −0.2−0.10.00.10.20.3 M0 M1 M2 M3b2 0.200.250.300.35 M0 M1 M2 M3c1 −1.6−1.5−1.4 M0 M1 M2 M3c2 −0.0035−0.0030−0.0025−0.0020 M0 M1 M2 M3c3 −0.50−0.45−0.40−0.35 M0 M1 M2 M3k 0.050.100.150.20 M0 M1 M2 M3f1 0.00.1 M0 M1 M2 M3f2Figure 6: Estimated coefficients for different models. Dot is median, uncertainty is 2.5% and 97.5% quantiles of posterior distribution. 5 Results for Real Data In this section, the model is applied to the observed data, with PGA as the target variable. Different models are used: •Model 0: Standard model, no spatial correlation, event and station terms. •Model 1: Similar to standard model, include spatial correlation of within-event/within- site residuals δWS •Model 2: no site terms, include spatial correlation of within-event residuals δW •Model 3: Like Model 1, but include cell-specific attenuation as in Kuehn et al. (2019) and Kuehn (2023). In addition, the spatial correlation parameters are also estimated for Models 1 and 2 from the residuals (i.e. using δWS for Model 1 and δWfor Model 2). Figure 7 shows the posterior distributions for the range and standard deviation of the spatial correlation structure. The spatial range differs between all models, with Model 2 (no site terms) 10 0.020.040.06 100 150 spatial range (km)posterior densityModel 1 Model 2 Model 3 full residual 255075 0.15 0.20 0.25 0.30 φcposterior densityModel 1 Model 2 Model 3 full residualFigure 7: Posterior distributions for practical spatial range and standard deviation ϕcwith real data. delta W range = 112.64km 0.000.030.060.090.12 0 100 200 300 400 h (km)γ(h) number of points 0 10000 20000 delta WS range = 83.82km 0.000.010.020.030.04 0 100 200 300 400 h (km)γ(h) number of points 0 10000 20000 Figure 8: Variograms for δWandδWS , aggregated for all events. having the largest value. Including cell-specific attenuation leads to a decreasing value of the spatial range; this makes sense, since spatial correlation of δWS has been interpreted as a proxy for path effects ( Anderson and Uchiyama , 2011; Kuehn and Abrahamson , 2020), and the cell-specific attenuation model already takes care of path effects. When estimated on residuals, the values are lower, similar to the simulation example, which is due due to ignoring uncertainty in event and site terms. One thing to note is that the spatial ranges for all models are quite long (on the order of about 100km), compared to other models (see Sgobba et al. (2023) for a comparison of effective ranges based on Italian data with other models; note that in Sgobba et al. (2023), the effective range is the distance where the correlation reaches a value of 0.05). The effective range estimated in Sgobba et al. (2023) is about 33km for PGA (see their Figure 7), which gives a practical spatial range of 21.71km, considerably shorter than estimated in this work. The reason is that during the regression, the spatial correlation structure is estimated for all records in the data set, aggregating over all events. On the other hand, dedicated spatial correlation studies typically select only well-recorded events, and estimate spatial correlation parameters separately for each event. The estimated length-scales/ranges often show large variation across events. Figure 8 shows variograms for residuals δWandδWS , aggregated for all events in the data set. The estimated spatial ranges based on these variograms (using the Mat´ ern model with ν= 1) are 83.82km for δWS and 112.64km for δW, which is in the same ballpark with the results from Figure 7. Figure 9 shows the posterior distributions of τ,ϕS2S, and of the standard deviation of the remaining residuals. As already seen in the simulation example, accounting for spatial correlation of residuals leads to a decrease in the value of τcompared to the base model. The ratio of variances is τ2 M0/τ2 M1= 1.82 (computed using the mean of the posterior distribution), which is similar to the values for the longer ranges in Figure 3. The models with spatial 11 102030 0.05 0.10 0.15τposterior densityModel 0 Model 1 Model 2 Model 3 204060 0.18 0.20 0.22 0.24 0.26 φS2Sposterior densityModel 0 Model 1 Model 3 50100150 0.10 0.15 0.20 0.25 0.30 φ0posterior densityfull residual Model 0 Model 1 Model 2 Model 3Figure 9: Posterior distributions for standard deviations with real data. 51015 0.15 0.20 0.25 0.30 0.35 0.40σcellposterior densityfull no corr Figure 10: Posterior distributions for standard deviations of cell-specific attenuation. correlation lead to lower station-to-station variability ϕS2S, which might indicate some trade- off between path and site effects. The standard deviations of the remaining residuals are quite different between the models, since they represent different quantities, depending on which random effects are included in the model. Figure 10 shows the standard deviation of the cell-specific attenuation coefficients, calculated both for a model including spatial correlation of δWS and without. Including spatial correlation leads to a decrease in σcell, meaning less variability in the cell-specific attenuation coefficients. This again indicates that both parts of the model target a path-specific effect and trade off in the model fit. Table 3 shows a summary of the standard deviations and some model selection criteria of the different models. In addition to the previously presented models, I also show results for a standard model (no correlation) without site effects, and a model including cell-specific attenuation, but no correlation. The second column in Table 1 shows which random effects are Table 3: Summary of results for models with different random effect structures. Shown are total and single-station standard deviations, and model selection criteria. See text for details. For WAIC and RMSE a lower value is better, for LL a higher value is better. Model RE σT σSS WAIC RMSE(LEO-CV) LL(LEO-CV) M0 δB,δS 0.341 0.248 -769.9 0.238 20.9 M0b δB 0.354 NA 918.1 0.332 -671.1 M1 δB,δS,dC 0.337 0.271 -2689.7 0.265 -202.9 M2 δB,δC 0.369 NA -445.4 0.337 -718.9 M3 δB,δS,δC, cell 0.309 0.243 -2716.1 0.238 9.4 M3b δB,δS, cell 0.304 0.230 -987.1 0.230 98.5 12 included in the model. The total standard deviation σTand the single-site standard deviation σSSare generally similar between models, with Model 1 showing slightly larger single-station sigma due to larger ϕSS. The models including cell-specific attenuation show smaller total sigma, since part of the ground-motion variability is accounted for by the variability in the attenuation effects. Table 3 also shows model comparisons based on criteria that evaluate how well the models can generalize to new, unseen data. Shown are the values of the widely applicable information criterion (WAIC, Watanabe , 2010, 2013), which is an information criterion for Bayesian models. Similar to the Akaike information criterion (AIC, Akaike , 1998), WAIC is a method to assess the out-of-sample prediction error of a model fit ( Vehtari et al. , 2017), and like AIC it penalizes model complexity. The models including spatial correlation clearly outperform their standard counterparts in terms of WAIC, with the cell-specific model leading to a small improvement. However, one should note that the model with correlation but excluding systematic site terms has larger WAIC than the standard model (no correlation but including event and site terms). Thus, adding site terms provides a much larger improvement in predictive capability than adding a spatial correlation structure. These results are not really surprising, since models with more random effects essentially use more information for prediction of new data points. Ming et al. (2019) provide a very good overview and discussion about the benefits of adding spatial correlation models for within-event prediction. WAIC is a measure of predictive accuracy, conditional on all random effects. Hence, it is akin to leaving out one record at a time and estimating its predictive distribution, including event terms, station terms, and other terms from the fit on the remaining records. However, the main use of GMM predictions is for PSHA, where one has to predict for new events, i.e. one cannot leverage observed records together with the spatial correlation model to improve predictions. Hence, Table 3 also shows the results of an experiment where all records from one event are left out as testing data, the different models are estimated based on the remaining training data, and the test data is compared to model predictions. This procedure is called ‘leave earthquake out-cross validation’ (LEO-CV). Predictions still include random effects from the training data, which in this case mean site terms and cell-specific attenuation coefficients, if applicable to the model. Both the root-mean-square error (RMSE) of the test residuals, as well as the total log-likelihood (LL), which is the sum of the LL values of all test data points, are calculated. The LL of each test data point is calculated as the value of the probability density function of the predictive distribution, which includes aleatory variability, but also uncertainty in the coefficients, hyperparameters and random effects. For LEO-CV, the ranking of the models is different, and the models including spatial correlation have worse predictive performance than their corresponding base models, both in terms of RMSE and LOO. Both WAIC and LEO-CV are calculated for events with at least 70 recordings. As shown in Figure 7, there is uncertainty associated with the spatial range parameters, quantified by their posterior distributions. The effect on the correlation is shown in Figure 11, as uncertainty bands corresponding to the 2.5% and 97.5% quantiles of the posterior distribution of the spatial range. In general, accounting for uncertainty can be important, and it should be investigated whether the results from a forward application of a spatial correlation model are sensitive to uncertainty in its parameters. However, given the large between-model variability in spatial correlation the within-model uncertainty is probably of minor importance. 13 0.000.250.500.751.00 0 50 100 150 200 distance (km)correlationFigure 11: Posterior distribution for spatial correlation. 5.1 Non-stationary Correlation Model Sgobba et al. (2023) argue that a simple stationary correlation model is not appropriate due to residual trends at longer distances. They achieve a non-stationary model by detrending the residuals; the resulting correlation model has a smaller spatial range. When including spatial correlation in the GMM regression, one cannot detrend the residuals, but needs to incorporate the non-stationarity directly in the correlation model. While it is generally well recognized that correlation of ground motion is non-stationary ( Chen and Baker , 2019; Chen et al. , 2021; Infantino et al. , 2021; Kuehn and Abrahamson , 2020), direct modeling is rare ( Bodenmann et al. , 2023; Kuehn and Abrahamson , 2020; Liu et al. , 2022). One can include covariates in the SPDE model ( Ingebrigtsen et al. , 2014), which allows to use a model similar to the one proposed by Kuehn and Abrahamson (2020), where the spatial range depends on the source- to-site distance. Thus, I now investigate a model where the spatial range is modeled as ℓ=( exph a1+a2Repi Rrefi Repi≤Rref exp [a1+a2] Repi> R ref(13) The reference distance Rrefis introduced to prevent the spatial range to reach ever larger values with increasing source-to-site distance. Here, its value is set to Rref= 80km, based on Sgobba et al. (2023). The implementation in INLA follows Chapter 5 of Krainski et al. (2019). Figure 12 shows the estimated spatial range dependent on epicentral distance, together with the estimate of the stationary model (Model 1 in previous section). Now, the spatial range close to the source is small and increases to a large value, much longer than for the stationary model. The effect is shown in Figure 13, where the correlation around an observation is displayed, color-coded by the value of the correlation. For an observation close to the source, the correlation drops rapidly, and most points are basically uncorrelated. By contrast, more points are correlated far away from the source. One can also see that the shape of the correlation is not symmetrical, but skewed away from the source. The stationary model is of circular shape and is the same everywhere, regardless of distance to the source. The WAIC value of the non-stationary model is -2984.2, which is lower than for the station- ary model (cf. Table 3). This indicates that the non-stationary model improves within-event prediction. 14 050100150200250 0 50 100 150 epicentral distance (km)spatial range (km)nonstationary stationaryFigure 12: Spatial range for non-stationary model. Stationary −200−1000100200 −200 −100 0100 200 xy 0.000.250.500.751.00 Non−stationary − Site 1 −200−1000100200 −200 −100 0100 200 xy Non−stationary − Site 2 −200−1000100200 −200 −100 0100 200 xy Figure 13: Correlation in space for stationary and non-stationary models. The range for the non-stationary model depends on epicentral distance, so is smaller for case one (station has location ⃗t= [10 ,0]T) than for case two (station has location ⃗t= [100 ,0]T). For the stationary model, the correlation is independent of station location. The epicenter has location ⃗t= [0,0]T. 15 6 Discussion This paper shows how to estimate empirical GMMs including spatial correlation structure of residuals in a Bayesian way using INLA. In the model, the spatial effects are added as an ad- ditional random effect which is estimated during regression, together with all other parameters and hyperparameters. The main question this paper raises is whether one should actually incorporate spatial correlations during a GMM regression. The main use of a GMM is prediction for PSHA, which requires calculation of the predictive distribution for records from unobserved events. As seen in Table 3, the model including spatial correlation performs worse in terms of RMSE and LL when predicting for observations from left out events. On the other hand, from a physical/theoretical point of view the residuals should be spatially correlated, so it makes sense to incorporate this aspect in the model. In the end, while prediction is important, it is also important to make the model as realistic as possible. It is reassuring that the base model without spatial correlation performs well from a predictive perspective, as most published GMMs are of this kind. If one considers to incorporate spatial correlations into GMM development, one should investigate whether the inclusion provides a benefit in terms of model performance, as well as produces sensible results. One advantage of estimating the parameters of the spatial correlation structure during the regression is that the uncertainties of all parameters (coefficients, latent effects and hy- perparameters) is taken into account. In the simulation example, the values of spatial range and associated standard deviation are biased when estimated from residuals. Typical spatial correlation models are estimated from within-event residuals δW, in which case the range is conceptually a mix between the range of δWS andδW, so one would not necessarily expect that the range is the same. The regression model including spatial correlations leads to estimates of the spatial range which are longer than those from published models. This might be due to the fact that the regression model takes into account all records from all events, and provides a spatial range based on the accumulated data (cf. Figure 8). Spatial correlation models are typically based on a smaller subset of well-recorded events, which might lead to different results. One potential solution could be to estimate a different spatial range parameter for each event during the regression, but this might lead to problems for events that do not have many recordings. In addition, in a forward application one would need to take the variability between events into account; one can use the (weighted) mean, but the spread of different range values should not be neglected. Another possibility for the longer range is that it compensates for trends in the residuals due to unmodeled effects, such as nonergodic path effects. The inclusion of a cell-specific attenuation model, which is designed to account for path effects, leads to a reduction of the spatial range by about 20km (cf. Figure 7). This effect can also be directly incorporated into the correlation structure via a non-stationary correlation function. A model based on Kuehn and Abrahamson (2020) and Liu et al. (2022), where the correlation depends on source-to-site distance, leads to smaller estimates of the spatial range close to the source, with increasing range for longer source-to-site distances. This model was proposed to account for path effects, and leads to a improved predictive performance. Alternatively, spatial correlation models for path effects have been proposed that are based on angular distance between source-to-site azimuths (Bodenmann et al. , 2023; Walling and Abrahamson , 2012). 16 While spatial correlation models are generally estimated from within-event residuals δW, empirical GMMs often partition residuals into within-event/within-site residuals δWS and site terms δS, and possibly additional terms ( Al-Atik et al. , 2010). In general, inclusion of site terms provides a strong benefit in reducing aleatory variability and improving predictive performance (cf. Table 3). Accounting for site terms should be done during the regression stage ( Sahakian et al., 2018), so incorporating spatial correlations should be done in a model including site terms (again, cf. Table 3). This needs to be kept in mind when assessing and using such a model (Stafford et al. , 2019). For example, if one uses a spatial correlation model that is estimated for within-event/within-site residuals δWS , e.g. to sample ground-motion maps for a regional risk analysis, then the sampled values do not account for spatial correlation of site terms δS. This correlation should be included in a separate step, using a spatial correlation model for site terms (e.g. Chao et al. , 2020; Parker and Baltay , 2022). The regression model with spatial correlation includes a noise term δWS 0, which is uncor- related between records, and in geostatistical parlance is called a nugget effect. In general, one would expect a continuous ground-motion field, i.e. the spatially uncorrelated noise term should be zero. However, there is probably small-scale variability which cannot be resolved by the spa- tial station distribution. There can also be measurement error due to different instruments, which can contribute to unexplained variability. The model as implemented in this work is estimated via Bayesian inference using INLA. It should be seen as a possible alternative to the algorithm presented by ( Ming et al. , 2019) (not a replacement) for analysts who prefer to use Bayesian inference. INLA can be used for models with a linear predictor, which is somewhat restrictive, though often the basic functional form of a GMM can be written as a linear model, like the ITACA18 model used as the basis in this work. Ming et al. (2019) provide an approximation to make their scoring algorithm work for non-linear models. For the INLA model, package Inlabru ( Bachl et al. , 2019) provides an iterative method to incorporate nonlinearities, which could be used for more complicated models. INLA has been used to implement nonergodic GMMs based on spatially varying coefficient models ( Kuehn , 2023; Lavrentiadis et al. , 2022). Thus, the regression models including spatial correlation of residuals can be extended to include nonergodic terms such as systematic, spa- tially varying event and site terms. In this work, I include systematic path effects based on cell-specific attenuation models, which changes both the spatial correlation structure (practical spatial range and standard deviation) and the parameters of the cell-specific model. Including spatially varying event and site terms is unlikely to have a large effect on spatial correlations of residuals, but it can be advantageous to have a full model to account for all uncertainties and trade-offs in a model at the same time. The SPDE approach implemented for spatial fields in INLA fixes the smoothness parameter νin the Mat´ ern correlation function (Equation 9) to ν= 1 for two-dimensional fields ( Lindgren and Rue , 2015). One could use a different value for ν, which would require some approxima- tions ( Lindgren et al. , 2011), but be more in line with most models of ground-motion spatial correlation,which typically use a exponential correlation function ( ν= 0.5). Alternatively, one can treat νas a parameter to be estimated ( Bolin and Kirchner , 2020; Xiong et al. , 2022). Differences between spatial fields using the Mat´ ern with ν= 1 and ν= 0.5 are likely small (Kuehn , 2021a), but it might be worth investigating. 17 Code and Data Availablity R-code to fit regression models including spatial correlations, as well as to generate some of the figures in this paper, can be found at https://github.com/nikuehn/GM_Spatial_ Correlation , archived at https://doi.org/10.5281/zenodo.11070746 (Kuehn , 2024). The flatfile for the data used in the ITA18 model can be found at https://shake.mi.ingv.it/ ita18-flatfile/ (Lanzano et al. , 2022). Acknowledgements The work is carried out in the computer environment R(R Core Team , 2021), using package R-INLA version 22.12.16. Plots are made with package ggplot2 (Wickham , 2016). Other packages used are Inlabru (Bachl et al. , 2019), tidyverse (Wickham et al. , 2019), gstat (Gr¨ aler et al. , 2016; Pebesma , 2004), and sf(Pebesma , 2018). I would like to thank Melanie Walling for multiple discussions about spatial correlation models with INLA. Some general ideas in this work came out of the nonergodic working group at UC Berkeley/UCLA. Partial support of Pacific Gas & Electric Company and California Department of Transportation are gratefully acknowledged. Any opinions, findings, and con- clusions or recommendations expressed in this material are those of the author and do not necessarily reflect those of the sponsors.