A nonrgodic ground-motion model explicitly takes systematic local source, path and site e ects on the predicted ground-motion into account. With an increasing number of ground-motion records, it s possible to estimate these e ects. Landwehr et al. (2016) proposed a varying coecient model as a tool to estimate nonergodic ground-motion models, based on Gaussian processes. Gaussian processes are computationally expen- sive, so for large data sets, some approximations have to be made. Here, we compare di erent Bayesian implementations of varying coecient models, using the probabilistic programming language Stan (Carpenter et al., 2017), and the integrated nested Laplace approximation (INLA) (Rue et al., 2009). The models are used to t nonergodic models on the California subset of the NGA-West2 data set (Ancheta et al., 2014). We nd that both implementations lead to very similar results, both in the estimated parameters and the predicted ground motions. 1 Introduction The goal of probabilistic seismic hazard analysis (PSHA) is to estimate the expected annual rate of exceedance of a ground-motion parameter of interest for a particular site. For that purpose, in PSHA one combines a seismic source model (SSC), which describes the temporal and spatial occurrence of events, with a ground-motion model (GMM), which models the distribution of the ground-motion parameter, given a particular earthquake scenario. The value of the standard deviation of aleatory variability associated with the ground-motion distribution controls the slope of the hazard curve, in particular at low exceedance rates (Bommer and Abrahamson, 2006). GMMs are typically estimated empirically using a regression approach. Since empirical data sets are sparse, one often relies on the ergodic assumption (Anderson and Brune, 1999) in some way, which states that the ground-motion distribution over time at one station is the same as the ground-motion distribution averaged over many stations in a particular region. Due to the increasing number of ground-motion recordings over the last decades, it is possible to relax the ergodic assumption, for example by estimating systematic station terms (e.g. Atkinson, 2006; Rodriguez-Marek et al., 2014). In addition, recent GMMs often include regional adjustment ∗kuehn@ucla.edu 1 terms (e.g. Sta ord, 2014; Kotha et al., 2020; Sedaghati and Pezeshk, 2017; Parker et al., 2020; Kuehn et al., 2020) to account for regional variations in site and distance scaling. These adjustments, however, are still applicable on a large regional scale, and thus employ the ergodic assumption over the regions underlying the models. Di erent methods have been proposed to further relax the ergodic assumption with respect to path e ects (Sung and Lee, 2016; Anderson and Uchiyama, 2011; Kuehn and Abrahamson, 2020). Dawood and Rodriguez-Marek (2013) proposed to model source-to-site speci c path e ects as the sum of atteuation over small cells, which was extended by Kuehn et al. (2019) using Bayesian inference to account for the uncertainties associated with the cell-speci c attenuation coecients. Landwehr et al. (2016) proposed a fully nonergodic GMM for California, based on a varying coecient model (VCM) (Gelfand et al., 2003; Bussas et al., 2015, 2017), where the coecients of the GMM are functions of source and site location. In the VCM of Landwehr et al. (2016), however, the anelastic attenuation parameter was modeled as dependent only on the event location, and thus did not fully capture path e ects, which should depend on both event and station locations. Abrahamson et al. (2019) combined the models of Landwehr et al. (2016) and Kuehn et al. (2019) to calculate nonergodic PSHA for three sites in California. They found that for sites with empirical data in the vicinity, di erences between the (mean) ergodic and nonergodic hazard curves can be quite large, while the mean nonergodic hazard curve for a site with no empirical data close-by is similar to the mean ergodic hazard curve, but with much wider fractiles, due to the epistemic uncertainty associated with the nonergodic terms. Any relaxation of the ergodic assumption means that some pars of the GMM are estimated with few data points; for example, when one calculates a regional adjustment to a coecient, the region under study has less data than the overall data set, and the coecient is thus less well constrained. It is thus very important to accurately estimate the epistemic uncertainties associated with any nonergodic term, since otherwise the hazard would be underestimated (Abrahamson et al., 2019). In the VCM of Landwehr et al. (2016), the spatial dependence of the varying coecients was modeled using Gaussian processes (GPs) (Rasmussen and Williams, 2006). GPs are a attractive choice for spatially varying parameters since they provide a exible framework to model the spatial dependence, while also providing good quanti cation of the associated uncertainty. In Bayesian nonparametrics (Hjort et al., 2010), a GP can be thought of as a prior over functions. Hence, we do not need to know the exact functional form of the spatial dependence; it is estimated from data, with some constraints set by the form of the GP. The GP formulation also lends itself to a good quanti cation of uncertainty in ranges where data is sparse (Landwehr et al., 2016) { this uncertainty can be calculated as the posterior distribution of the function describing the spatial dependence. While GPs provide an elegant mathematical framework for estimating nonergodic GMMs and assess the associated uncertainty, they do, however, have the drawback that they do not scale well to large data sets. In particular, estimating the parameters of a GP (as well as predicting new data points using a GP) requires the calculation of the inverse of large covari- ance matrices, which is computationally expensive. For large data sets, this requires some approximations to the covariance matrix. In the future, the size of strong-motion data sets will continue to increase, as more data becomes available. This o ers many modeling opportunities, since more data allows one to better constrain nonergodic e ects. On the other hand, this also poses computational challenges, as models with spatial components are harder to estimate. It is therefore necessary to explore how nonergodic models can be t eciently for large data sets. Landwehr et al. (2016) estimated the parameters of their nonergodic GMM by maximizing the marginal likelihood of the data (cf. Chapter 5.4.1 in Rasmussen and Williams (2006)), which gives point estimates of the parameters. However, parameters are often associated with considerable uncertainty, which should not be ignored. Hence, instead of point estimates it is 2 Figure 1: Caption better to work with the posterior distribution of all parameters, and estimate the models via full Bayesian inference. Here, we compare di erent implementations of Bayesian nonergodic models (a combination of a VCM with the cell-speci c attenuation model of (Kuehn et al., 2019)). We estimate nonergodic models using three di erent methods: (1) a full Bayesian approach without approximation to the GP, estimated via Hamiltonian Monte Carlo (HMC) sampling (Neal, 2011; Betancourt, 2017; Betancourt and Girolami, 2015); (2) a full Bayesian approach using the nearest neighbour GP, (Datta et al., 2016a), which is an approximation too a full GP, estimated via HMC; and (3) a full Bayesian approach using the integrated nested Laplace approximation (INLA) (Rue et al., 2009). We compare the resulting nonergodic models, in terms of their estimated parameters, and their predictions. 2 Data The parameters of the nonergodic adjustment models are estimated using the Western US data from the model of Abrahamson et al. (2014). The dataset comprises 12018 records from 274 earthquakes and 1479 stations. For simplicity, we use the same data set for each of the three ergodic base GMMs (ASK14, BSSA14, CB14). The event and station locations are is shown in Figure 1. We estimate the nonergodic adjustment models for a pseudo-spectral acceleration (PSA) with a period of T= 0:02s. The median predictions for the ergodic base GMMs are calculated using the python package pygmm (Kottke, 2016). For the cell-speci c attenuation, we use the same cells as Kuehn et al. (2019), who used cells of about 2525km. The cell-speci c distances are calculated as a straight line from the closest point on the rupture to the station. Hence, the The coordinates of the events are the latitude and longitude of the hypocenter coordinates in degrees. 3 3 Nonergodic GMMs and PSHA The general approach to performing nonergodic PSHA is outlined in e.g. Abrahamson et al. (2019); Villani and Abrahamson (2015); Kuehn and Abrahamson (2020). The key di erence between an nonergodic and an ergodic ground-motion model is that the nonergodic GMM explicitly takes into account repeatable, systematic source, path and site e ects. These sys- tematic e ects are location speci c; hence, the median prediction of a GMM becomes location dependent nonergodic (~ x;~te;~ts) =ergodic (~ x) +feq(~te) +fpath(~te;~ts) +fstat(~ts) (1) wherefeq(~te);fpath(~te;~ts);fstat(~ts) are the location-speci c source, path and site e ects, respec- tively; and~teand~tsare the event and site coordinates. ergodic (~ x) is the prediction of an ergodic GMM, and ~ xis a vector of predictor variables such as magnitude, distance, time averaged shear wave velocity in the upper 30m VS30, and so on. For the current study, location-speci c nonergodic adjustment terms are calculated for the empirical GMMs of Abrahamson et al. (2014); Boore et al. (2014); Campbell and Bozorgnia (2014), hereafter called ASK14, BSSA14, and CB14. The nonergodic adjustment terms are based on a varying coecient model Landwehr et al. (2016); Bussas et al. (2017), and the cell- speci c attenuation model Dawood and Rodriguez-Marek (2013); Kuehn et al. (2019). We do not calculate adjustment terms for the model of Chiou and Youngs (2014), because it models the anelastic attenuation term as magnitude dependent, which is dicult to incorporate into the cell-speci c attenuation model. The nonergodic adjustment models are estimated via regression from California data. We use two main methods (di erent programs) to estimate the nonergodic models: (1) Bayesian inference via Markov Chain Monte Carlo (MCMC) sampling using the program Stan (Car- penter et al., 2017); (2) Bayesian inference vie the R-INLA package (Rue et al., 2009, 2017; Bivand et al., 2015). For Stan, we describe an exact GP model model, as well as the nearest neighbour Gaussian process approximation (Datta et al., 2016a). There are some di erences in the methods, described below. From a practical point of view, the main di erences are in the runtime of di erent approaches. The exact version of Stan takes about 120 hours for the California dataset, the nearest neighbour approximation runs for about 20 hours, and INLA runs about 2 hours. 3.1 Nonergodic Adjustment Models This section describes the formulation of the nonergodic adjustment models. We assume that we have an ergodic GMM, whose median prediction can be calculated by GMM =fGMM (~ c;~ x) (2) where~ xis a vector of predictor variables such as magnitude, distance and so on; and ~ cis a vector of coecients of the model. We want to make this GMM nonergodic, which means that we need to adjust the median prediction with systematic, location-speci c nonergodic adjustment terms. The adjustment terms are an earthquake-speci c constant; a site-speci c constant; and a cell-speci c anelastic attenuation term. Following Landwehr et al. (2016), the constants are modeled as spatially correlated (i.e., they are distributed according to a Gaussian process). The geometrical spreading is modeled as ergodic (i.e. not spatially varying); this is in line with most partially nonergodic models, which only regionalize the constant, linear site scaling, and anelastic attenuation. The cell-speci c anelastic attenuation coecient is modeled as in Kuehn et al. (2019), following Dawood and Rodriguez-Marek (2013). We also include an intercept, to account for any model misspeci cation. 4 Hence, the full model can be written as lnY=c0+fGMM (~ c;~ x) +ceq(~te) +cstat(~ts) fattn;GMM (RRUP)NCX j=1cca;jRj +B 0+W 0+S2S0 (3) Theceq=stat -terms are the location-speci c adjustment terms for evens and stations, and ~te=sare the event and station locations. The term fattn;GMM (RRUP) is the anelastic attenuation function in the GMM; it is taken out and replaced by the cell-speci c attenuation,PNC i=1cca;iRi, where cca;iis the coecient for cell i, andRiis the length of the path within the ith cell. We constrain the cell-speci c attenuation coecients to be positive, so the full attenuation term has a negative sign.B 0;S2S0;W 0are the remaining, spatially uncorrelated between-event, site-to-site, and within-event/within-site residual. These residuals are assumed to follow normal distributions with between-event variability 0, site-to-site variability S2S;0, and within-event/within-site variability 0. B 0N(0; 0) (4) S2S0N(0; S2S;0) (5) W 0N(0; 0) (6) The location-speci c constants ( ceq(~te) andcsite(~ts)) are distributed according to a Gaussian process prior, as in Landwehr et al. (2016) cGP(0;k(~t;~t0) (7) wherek(~t;~t0) is the covariance function that determines how the adjustment terms can vary. Since thecare adjustments to an ergodic base model, a zero mean function is appropriate. De- pending on the method to estimate the ergodic adjustment terms, we either use an exponential covariance function (for estimation with Stan) or a Mat ern covariance function (for estimation with INLA). In CB14, the anelastic attenuation coecient is applied to records wih RRUP>80km. We apply the same constraint to the cell-speci c attenuation model, if the base model is CB14. Hence, in this case the cell-speci c distance  Rjis zero if the total distance up that point is less than 80km. A quick word on notation and terminology: In the following, the adjustment terms for eventsceq, and stations cstat, are treated mathematically the same. Hence, we often refer to them in the text as c, which implies that the corresponding explanation applies to both event and station adjustments. Similarly we use ~tfor the locations when speaking generally about the model. Parameters that govern the distribution of other parameters (such as governs the distribution of event terms) are called hyperparameters. We refer to parameters whose distribution is controlled by hyperparameters as latent parameters or e ects. Latent parameters are for example the event and station terms, as well as the adjustment terms c. Below, we describe the derivations of the nonergodic models, using both Stan and INLA. We provide code to run the models using the computer environment R(R Core Team, 2020) at https://github.com/nikuehn/NonergodicGMMS_STAN_INLA . 3.2 Derivation of Nonergodic Parameters using Stan First, we describe how the parameters of the nonergodic adjustment model described in the previous section are estimated using full Bayesian inference, via MCMC sampling with the 5 program Stan (Carpenter et al., 2017). The covariance function we use is the exponential covariance function, as in Landwehr et al. (2016) k(~t;~t0) = 2exp j~t~t0j ` (8) where determines the range how much the adjustment terms can vary, and `is a spatial correlation length-scale. Depending on whether the GP is modeling the source or site-speci c adjustment terms, we use either the event or site locations in Equation (8). In general, one can integrate out all latent adjustment terms (non-ergodic adjustment terms as well as event and site terms) in Equation (3), which leads to a model where the target variable is distributed according to a multivariate normal distribution ~YN(;) (9) where i=c0+fGMM (~ c;~ xi)fattn;GMM (RRUP;i)NCX j=1cca;jRi;j (10) The covariance matrix is ij=ij 2 0+eq(i);eq(j) 2 0+stat (i);stat (j) 2 S2S;0+keq(~te;~t0 e) +kstat(~ts;~t0 s) (11) whereijequals 1 if i==jand zero otherwise, and eq(i) andstat(i) are indices connecting an event or station with record i. However, is of sizeNN, and is not sparse, so it is impractical to work with Equation (9), since one has to calculate the inverse of  multiple times during the estimation of the model parameters. If we do not integrate out the latent parameters, but treat them as parameters to be esti- mated, then the likelihood of the data factorizes, conditional on the latent e ects. Since the posterior distribution of the parameters ( ~, which encompasses all parameters and latent param- eters of the model) is proportional to the product of the likelihood and the prior distribution, we can write the full log-posterior distribution (up to a constant) as lnp(~jD) = lnp(c0; 0; S2S;0; SS;0; eq;`eq; stat;`stat;ca; ca) +NEX i=1lnp(B 0;ij0; ) +NSX j=1p(S2S0;jj0; S2S) + lnp(~ceqj~0;Keq) + lnp(~cstatj~0;Kstat) +NCX j=1lnp(cca;jjca;ca)NCX j=1ln(1cdf(0jca; ca)) +NX k=1lnp(ykjk+Beq(k)+S2Sstat (k)+ceq(k)+cstat (k); SS) (12) whereDdenotes the data, and p(jm;s) is the probability density function of a normal dis- tribution with mean mand standard deviation s(or a multivariate normal distribution if the argument is a vector). In Equation (12), the rst line describes the prior distribution for the interceptc0and the hyper-parameters of the model (i.e., parameters that govern the distribu- tion of the latent e ects). The second line describes the prior distributions for the systematic, non-spatially varying event and station terms. These are independent of each other, and are normally distributed with mean zero and standard deviation 0and S2S;0. 6 The third line in Equation (12) is the prior distribution of the systematic event and station terms, which is a multivariate normal distribution with mean zero and covariance matrix Keq=sat de ned by the covariance function (cf. Equation (8)). The covariance matrix Keqhas size NENE. The fourth line in Equation (12) is the prior distribution for the cell-speci c attenuation models. We follow the formulation of Kuehn et al. (2019), and assume that the cell-speci c attenuation coecients are distributed according to a truncated normal distribution. Here, we constrain the parameters to be positive, so the prior distribution is lower truncated at zero, and the second term is the required adjustment to the untruncated density. A Stan program consists of statements that increments the log-density for each component of the model, so Equation (12) is straightforward ti implement. For GPs, informative prior distributions are important, to constrain the range of the pa- rameters, and help to avoid possible trade-o s. For the intercept c0, we do not believe that it is very large, since the models should be approximately centered. Hence, we give it a normal prior with a small standard deviation c0N(0;0:1) (13) For the standard deviation parameters 0; 0; S2S;0, we use a lognormal distribution that discourages values below 0.2 and above 1.2. 0lognormal(ln(0 :5);0:5) (14) 0lognormal(ln(0 :5);0:5) (15) S2S;0lognormal(ln(0 :5);0:5) (16) The rationale for the prior distributions for the scale parameters of the GPs, , is to penalize too large values of , since these would be a strong deviation from the ergodic base model (Simpson et al., 2017). We do believe that the nonergodic values exist, but we need to see strong evidence in the data to model them. Hence, we use an exponential distribution as the prior for these parameters eqexponential(20) (17) statexponential(20) (18) Kuehn et al. (2020) compared prior distributions for regional partially nonergodic parameters based on an exponential distribution and a half-cauchy distribution, which is sometimes used as a prior for standard deviations in Bayesian hierarchical models (Gelman, 2006). They found that the exponential distribution gave better results. For the length-scale `, we use an inverse Gamma distribution as a prior distribution1, since this discourages both very small and very large length-scales. We scale the parameters of the inverse Gamma distribution such that length scales below 0.06 and above 1.2 are discouraged. This leads to the following prior distributions `eqinvgamma(3;0:5) (19) `statinvgamma(3;0:5) (20) For the cell-speci c attenuation coecients, we constrain them to be positive, and then use a minus sign in Equation (3). The cell-speci c coecients ccaare distributed according to a normal distribution with mean  caand standard deviation cca. The prior distribution for the mean attenuation is a lognormal distribution (which is why we constrain the parameters to 1based on https://betanalpha.github.io/assets/case_studies/gp_part3/part3.html 7 be positive) with median ln 0 :011, while the prior distribution for the standard deviation is a truncated Student-t distribution with degrees of freedom = 3. ccaN(cca; cca)T(0;) (21) ccaLN(4:5;0:5) (22) ccastudent-t(3;0;0:01)T(0;) (23) 3.2.1 Nearest Neighbour Gaussian Process In the Stan code of the model, the computational bottleneck is the calculation of the covariance matrices KstatandKeqand their Cholesky decompositions. The size of these matrices (in particular Kstat) leads to long runtimes of the Stan program. One way to overcome this issue is to use an approximation to the covariance matrices. Here, we use a nearest neighbour Gaussian process (NNGP) (Datta et al., 2016a,b) to approximate the original GP. The basic idea behind the NNGP is to represent the multivariate variables ~cas a directed acyclic graph (DAG) p(~c)p(c1)NY i=2p(cijPa[i]) (24) whereci=c(~ti), and Pa[i] is the set of parents for ci. In the NNGP, the size of Pa[ i] is shrunk to a small number (corresponding to the nearest neighbors of ~ti), leading to a sparse covariance matrix, which can be exploited computationally. For details on the NNGP approximation, see (Datta et al., 2016a,b; Finley et al., 2018; Zhang et al., 2018; Zhang, 2018). The calculation of the likelihood of the NNGP is taken from Zhang (2018). For our models, we use a number of 10 nearest neighbours for events and stations. We use the same prior distributions as in the Stan program without the NNGP approximation. 3.2.2 Running the Stan Programs We run the Stan programs using the command line interface cmdstan . We run four chains for each case, with 1000 warm-up samples and 1000 samples post warm-up. 3.3 Derivation of Nonergodic Parameters using INLA In the previous section, we described how the nonergodic parameters are calculated using Stan. Via Stan, the parameters are estimated using MCMC sampling, which means that the outcome of the program are samples of the full posterior distribution of all parameters. From these samples, one can assess summary statistics such as mean, median, standard deviation for each parameter, or correlations between di erent parameters. By contrast, INLA (Rue et al., 2009) provides a deterministic method to perform Bayesian inference by approximating posterior marginal distributions. The INLA method is implemented in the R-package R-INLA ( www.r-inla.org ) (Rue et al., 2017). For an introduction to INLA and example models, see e.g. Krainski et al. (2019); G omez-Rubio (2020). In particular, INLA can be used for eciently estimating spatial models combining latent Gaussian models with stochastic partial di erential equations (SPDEs) (Lindgren et al., 2011; Lindgren and Rue, 2015; Bakka et al., 2018). Using INLA for the estimation of spatial models requires a slight change of the covariance function. In Equation (8), we have used an exponential covariance function. The covariance function used in INLA is a Mat ern covariance function K(~t;~t0) = 22(1) ()(j~t~t0j)K(j~t~t0j) (25) 8 33363942 −125.0 −122.5 −120.0 −117.5 −115.0 Longitude (deg)Latitude (deg)Figure 2: Mesh used for the INLA models. Red dots are the stations, blue dots are the events. where  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). The default value in INLA for GPs with two-dimensional coordinates is = 1 (Lindgren and Rue, 2015; Bakka et al., 2018). Other values are possible, but require an additional approximation. Spatial modeling with INLA approximates a spatial process (in our case, the systematic event and station terms) with a basis function approximation (Lindgren et al., 2011; Krainski et al., 2019). The basis functions are evaluated on nodes of a triangular mesh. We generate a mesh that covers the state of California and the observed events and stations, We use a mesh that encompasses the observed data points, with an outer bound of 0.5 degrees. The maximum triangle edge length is 0.05 degrees, and 0.25 degrees in the outer bound. The nal mesh consists of 17794 nodes and is shown in Figure 2. Based on this mesh, we t the model and estimate the parameters. The parameters are the same as before, but the systematic terms now are distributed according to a GP with Mat ern covariance with = 1. The prior speci cation is di erent compared to the Stan implementation. In INLA, the normal distribution is parameterized in terms of mean and precision, which is de ned asprec =1 2(i.e. precision is one over variance). For the variance components of the non-spatially varying terms ( 0; 0; S2S;0), we place Gamma distribution on the corresponding precision parameters. The prior distributions have shape parameter 2 and rate parameter 0.5 lnprec G(2;0:5) (26) lnprec 0G(2;0:5) (27) lnprec S2S;0G(2;0:5) (28) Thus, the prior distributions for the precisions have a mean of 4, which corresponds to a value of the standard deviation of 0.5, which is similar to the prior speci cation that was used in the Stan implementation. We again assume that the base models should be approximately centered, and assign a normal distribution with mean zero and precision 100 to the intercept c0. The prior values for the parameters controlling the spatially varying parameters are so-called penalized complexity (PC) priors (Simpson et al., 2017; Fuglstad et al., 2019; Franco-Villoria 9 Table 1: Model selection criteria for ergodic and nonergodic INLA models. GMM WAIC (ergodic) WAIC (VCM) CPO (ergodic) CPO (VCM) ASK14 20044.2 18982.8 10074.4 9541.43 BSSA14 19923.3 18857.9 10016. 9480. CB14 19922.4 19361.2 10014.2 9729.19 et al., 2018). These priors shrink the model towards a base model without a spatial component (i.e. they penalize the extra complexity of the spatial components). The simpler base model has zero variance and in nite range, where the range is de ned as `=p 8 . The PC prior can be speci ed through the tail probabilities P( > 0) = 1andP(`<` ) = 2. We assume the same values for both the spatial model associated with the event and station terms. We use the following prior speci cation for the PC-prior of the spatially varying coecients P( eq=stat>0:23) = 0:01 (29) P(`eq=stat<0:0:45) = 0:9 (30) These values are derived from the prior distributions on the spatial hyperparameters used in the Stan program. For details on penalized complexity priors in the context of VCMs, see Franco-Villoria et al. (2019). The anelastic attenuation term is implemented in the INLA model as a xed e ect (corre- sponding to cain Equation (22)) and cell-speci c adjustment coeents. The prior distribution for the xed e ect is a normal distribution with mean 0:01 and precision 10000. The cell- speci c adjustment coecients are implemented as generic random e ects, since they can be written as NCX i=1cca;iRi=Z~ cca (31) whereZis a matrix of size NNC, whosei;jth entry corresponds to the fraction of path in celljfor recordi. The R-INLA package contains a module to specify random e ects of the form Z~ c, where Zis the design matrix of the random e ects, and ~ cis a vector of random e ects. 4 Results In this section, we brie y summarize results of the di erent models. We compare the poste- rior distributions of the estimated parameters, for the di erent base models and the di erent estimation procedures. For the Stan based models, we assess the posterior distribution via the post warm-up samples from the MCMC chains, while the outcome of INLA are approximations of the marginal posterior distributions for the di erent parameters. For comparison, we also estimated a model without spatially varying coecients with INLA, i.e. a model that only ac- counts for an intercept, event terms Band station terms S. In traditional GMM terminology, such a model would be called partially ergodic, since it accounts for systematic station terms, but for brevity we call this model ergodic. In Table 1, we compare the model t of the ergodic INLA models and the INLA VCMs, based on the widely applicable information criterion (WAIC, Watanabe, 2013), and the Conditional Predictive Ordinate (CPO, Pettit, 1990). For both criteria, lower values indicate a better t. In Table (1), the VCMs perform better than the ergodic models, which is in line with previous results by e.g. Landwehr et al. (2016); Kuehn et al. (2019), and also aligns with our expectation that moving from an ergodic to nonergodic models should lead to model improvement. 10 -���-���-�������������������������� �������������� ������ �������� ���� ���� ���� �������Figure 3: Posterior density of the intercept of nonergodic models for di erent base GMMs and estimation procedures. For comparison, the intercepts of ergodic models are also shown. Overall, the results are quite similar between the full Stan model, the NNGP Stan model, and the INLA model. Depending on the parameter, there can be some di erences between the adjustment models based on the di erent base models. This is not really surprising, since the NGA West 2 models can be quite di erent for small magnitudes, and thus the event and station terms for individual events/stations can be strongly a ected by di erent scaling of the base models. Figure 3 shows the posterior distributions of the intercept c0. In general, all models show similar behavior. We can also see that the intercept is di erent between the ergodic and nonergodic models. This can be due to the fact that an ergodic models can be biased towards regions where data is more abundant. One can also see that the posterior distribution for the VCMs is a bit wider than for the ergodic model. This is not surprising, since the model is distributed over more terms, which leads to larger uncertainty. Figure 4 shows the posterior distributions of the standard deviations of the non-spatially varying terms ( B 0,S0,W 0), i.e. of 0, S2S;0, SS;0. Again, we also show the results of an ergodic INLA run for comparison. As expected, we get a reduction for the nonergodic models compared to the ergodic runs. The largest decrease occurs for S2S, which is line with results of Chao et al. (2020). The results between the di erent approaches and base models are quite similar, in particular for within-event/within-station standard deviation 0, and site-to-site standard deviation S2S;0. For CB14, there is less of a reduction in 0, which is probably due to the fact that the anelastic attenuation term only applies for distances RRUP>80km in CB14. For the between-event standard deviation 0, the di erences are a bit larger. On reason could be that the between-model variability of the three base models is larger a small magnitudes, where the bulk of the data is. Overall, though, the di erences between the Stan results and INLA results is quite small. In particular, the full Bayesian model and the NNGP model produce almost identical results. The posterior distributions of the hyperparameters describing the spatially varying event and station terms are shown in Figure 5 and 6. The standard deviations of the systematic event terms ( eq) are quite consistent, but there are some di erences in the posterior distributions of the length scales. In particular, the INLA models lead to larger length scales. This is probably due to the di erent prior speci cation, as the PC prior used in INLA tries to shrink the estimates towards a simpler base model with in nite spatial range. A similar picture emerges for the hyperparameters associated with the systematic station terms, with INLA models having larger length scales. In this case, the models based on CB14 also have lager systematic standard deviations stat. This might be due to the fact that the anelastic attenuation in CB14 only applies to records with distances R> 80km. For the other base models (ASK14 and BSSA14), 11 (�) �������������������������������������������� ϕ������������� ������ �������� ���� ���� ���� ������� (�) ������������������������������������ τ�������� (�) ���������������������������������������� ϕ������������Figure 4: Posterior density of the standard deviations of the non-spatially varying terms, for di erent base GMMs and estimation procedures: (a) aleatory within-event/within-site variabil- ity 0; (b) aleatory between-event variability 0; (c) epistemic aleatory between-site variability S2S;0. For comparison, the distributions for ergodic models are shown as well. (�) ��� ��� ��� ��� ������������ ψ�������������� ������ �������� ���� ���� (�) � � � � ������ ℓ��������� Figure 5: Posterior densities of (a) the standard deviation; and (b) the length scale associated with spatially varying event terms. 12 (�) ��� ��� ��� ��� ��� �������������� ψ���������������� ������ �������� ���� ���� (�) ��� ��� ��� ��� ��������� ℓ�����������Figure 6: Posterior densities of (a) the standard deviation; and (b) the length scale associated with spatially varying station terms. (�) ������������������������������������������������������������� ψ�������������� ������ �������� ���� ���� (�) ��������������������������������������������������������� μ���������� Figure 7: Posterior densities of (a) the standard deviation of the cell-speci c atenuation coef- cients; (b) the mean anelastic attenuation. Vertical bars in (b) are the anelastic attenuation coecients of the ergodic base models. some of those records are a ected by the cell-speci c attenuation, while in CB14 some of the systematic e ects might be mapped into station terms. The parameters controlling the anelastic attenuation are shown in Figure 7. We can see a similar picture as before, with almost identical results between the full Stan model and the NNGP approximation, and only small di erences on the INLA models. As already observed by Kuehn et al. (2019), the di erences for cell-speci c parameters based on ASK14 and BSSA14 are almost negligible. Similar to other parameters the CB14 model shows a slightly di erent behavior, due to its di erent anelastic attenuation model. Figure 7(b) shows the mean atten- uationca, and for comparison the values of the (absolute) ergodic attenuation coecients. These are di erent, due to the fact that the majority of the data is from Southern California, so the ergodic attenuation coecients are biased towards Southern California. Kuehn et al. (2019) found that the ergodic attenuation coecient is a weighted average of the cell-speci c attenuation coecients, with weights corresponding to the number of paths in each cell. Figure 8 shows the di erence in the cell-speci c attenuation coecients, based on the full Stan model, the NNGP model, and the INLA model. As previously seen, the full Stan model and the NNGP Stan model give very similar results. The INLA coecients generally follow a similar trend; however, we can see that there are a few cell-speci c coecients that are negative, which leads to a positive anelastic attenuation, which is unphysical. For the Stan 13 ●●● ●●●●●●●●●●●●●● ●●●● ●●●● ● ●●●● ● ●● ●● ●●● ●●●● ●●●● ●●● ● ●●●●● ● ● ●● ●● ●●● ●●●●●●●● ● ●● ●● ●●●●● ●● ●●● ●●● ●● ●● ●●● ● ●●●● ●●● ●●●● ●●●●●●● ●● ●● ●● ●●●●●● ●●● ●●●●●● ●●●●●●●●●● ●● ●● ● ●●● ●●● ●●●●●●● ●●● ●●● ●● ●●● ● ●●● ●●● ●●●●●●● ●●● ●●●●● ● ●● ● ●●● ●●●●● ●●● ●●● ●●●●● ● ●●●●●●●●●● ● ●●● ●●●● ●● ●● ●●● ●● ●●●● ●●● ●●●●●●● ●●●●● ● ●●●●●●● ●● ●● ●●● ●●●●● ●●●●●●● ●●● ●●● ● ● ●●●●●●●●●●● ●● ●● ●● ●● ●● ●●●●●●●●●●●● ●● ●● ●● ●●● ●●●●●●●●●● ●● ●●●●●● ●●●● ●● ●●●● ●●●●●●●●●● ●●●●●●●● ●●●●○○○ ○○○○○○○○○○ ○ ○○○ ○○○○ ○○○○ ○ ○○○○ ○ ○ ○ ○○ ○○○ ○○○○ ○○○○ ○○○ ○ ○ ○○○○ ○○ ○○ ○○ ○○○ ○○ ○○○○ ○○ ○ ○○ ○○ ○○○○○ ○○ ○○○ ○○○ ○○ ○○ ○ ○○ ○ ○○○○ ○○○ ○○○○ ○○○○○○○ ○○ ○○ ○○ ○○○○○○ ○○○ ○○ ○○○○ ○○○○○○○○ ○○ ○○ ○○○ ○○○ ○○○ ○○ ○○ ○○○ ○○○ ○○○ ○○ ○○○ ○ ○○○ ○○○ ○○○○○○○ ○○○ ○○○○○ ○ ○○ ○ ○○○ ○○○○○ ○○○ ○○○ ○○○○ ○ ○ ○○○ ○○○ ○○○○ ○ ○○○ ○○○○ ○○ ○○ ○○ ○ ○○ ○○○○ ○○○ ○○○○○○○ ○○○○○ ○ ○○ ○○○○○ ○○ ○○ ○○○ ○○○ ○○ ○○○ ○○○ ○ ○○○ ○○○ ○ ○ ○○○○○○○○○○○ ○○ ○○ ○○ ○○ ○○ ○○○○○○○○○○○○ ○○ ○○ ○○ ○○○ ○○○○○○ ○○○○ ○○ ○○○ ○○○ ○○○○ ○○ ○ ○ ○○ ○○○○○○○ ○○○ ○○○○○○○○ ○○○○ ���� ���� ���� �������������������� ���(���� ���� ����� )���●���� ○���� ●●● ●●●●●●●●●●● ●●● ●●●● ●●●● ● ●●●● ● ●● ●● ●●● ●●●● ●●●● ●●● ● ●●●●● ● ● ●● ●● ●●● ●●●●●●●● ●●● ●● ●●●●●●● ●●● ●●● ●● ●● ●●●● ●●●● ●●● ●●●● ●●●●●●● ●●●● ●● ●●●●●● ●●● ●● ●●●● ●●●●●●●● ●● ●● ●● ● ●●● ●●● ●●●●●●● ●●● ●●● ●● ●●● ● ●●● ●●● ●●●●●●● ●●● ●●●●● ● ●● ● ●●● ●●●●● ●●● ●●● ●●●● ● ● ●●●●●●●●●● ● ●●● ●●●● ●● ●● ●●● ●● ●●●● ●●● ●●●●●●● ●●●●● ● ●●●●●●● ●● ●● ●●● ●●●●● ●●●●●● ● ●●● ●●● ●● ●●●●●● ●●●●● ●● ●● ●● ●● ●● ●●●●●●●●●●●● ●● ●● ●● ●●● ●●●●●●●●●● ●● ●●●●●● ●●●● ●● ● ●●● ●●●●●●●●●● ●●●●●●●● ●●●●○○○ ○○○○○○○○○○ ○ ○○○ ○○○○ ○○○○ ○ ○○○○ ○ ○ ○ ○○ ○○○ ○○○○ ○○○○ ○○○ ○ ○ ○○○○ ○○ ○○ ○○ ○○○ ○○ ○○○○ ○○ ○ ○○ ○○ ○○○○○ ○○ ○○○ ○○○ ○○ ○○ ○○○○ ○○○○ ○○ ○ ○○○○ ○○○○○○○ ○○ ○○ ○○ ○○○○○○ ○○○ ○○ ○○ ○○ ○○○○○○○○ ○○ ○○ ○○○ ○○○ ○○○ ○○○○ ○○○ ○○○ ○○○ ○○ ○○○ ○○○○ ○○○ ○○○○○○○ ○○○ ○○○○○ ○ ○○ ○ ○○○ ○○○○○ ○○○ ○○○ ○○○○ ○ ○ ○○○○○○ ○○○○ ○ ○○○ ○○○○ ○○ ○○ ○○ ○ ○○ ○○○○ ○○○ ○○○○○○○ ○○○○○ ○ ○○ ○○○○○ ○○ ○○ ○○○ ○○○ ○○ ○○○○○○ ○○○○ ○○○ ○○ ○○○○○○○○○○○ ○○ ○○ ○○ ○○ ○○ ○○○○○ ○○○○○ ○○ ○○ ○○ ○○ ○○○ ○○○○○○ ○○○○ ○○ ○○○ ○○○ ○○○○ ○○ ○ ○ ○○ ○○○○○○○ ○○○ ○ ○○○○○○○ ○○○○ ���� ���� ���� �������������������� ���(���� ���� ����� )���●���� ○���� ●●●●●●●●●●● ● ●●● ●●●● ●●●● ●●●●●● ●●●●●●●●●● ●●●●●●●●● ●● ●● ● ●●●●●●●●● ●● ●●●●●●●●●● ● ●● ●●●●●●● ●●●●●●● ●●●●●● ●●●●● ●●●●● ●● ●●●●●●●●● ●●●●●●●●● ●● ●● ●●●●●●●●●●●● ●●●● ●●● ●● ●●●●●●● ●●●●● ●●●●●●●●●●●●● ● ●●●●● ●●●●●●●●●● ●● ●●●●● ●● ● ●●●●●●● ● ●●● ●●●● ●● ●●●● ●●●●●●●●●●●●●●●● ● ●● ●●●●●●●●●●●● ●● ●●●●●●●●●● ●●●●●●● ●●● ● ●●●●●●●●●●● ●●●● ●● ●●●●●●●●●●●● ● ●●● ●●●●●●●● ●●●● ●● ●●●●●●● ●●●●●●●●●●●●●● ●●●●○○○ ○ ○○○○ ○○○○ ○○○ ○ ○ ○○ ○○○○ ○○ ○ ○ ○○ ○ ○○○○ ○○ ○○○ ○○○○ ○ ○○○○ ○○ ○○ ○ ○○○○ ○○ ○○○ ○○ ○○○○○ ○○ ○○○ ○ ○○○○ ○○○○ ○ ○ ○○○ ○○ ○○○○○○ ○ ○○○○ ○ ○○ ○○ ○○ ○○ ○○○ ○○○○ ○○ ○○○○○○○ ○○ ○○ ○○ ○○ ○ ○○○○○○ ○ ○○○○ ○○○ ○○ ○○ ○○○ ○○ ○ ○○○○ ○○○ ○○○○ ○○○ ○○○ ○ ○○ ○○ ○ ○○ ○ ○○○○○○○ ○○ ○ ○○○○ ○○ ○ ○○ ○○○○○ ○ ○○○ ○○ ○○ ○○ ○○○ ○ ○○○○○○○○○ ○○○○○ ○○○ ○○ ○○ ○ ○○○○ ○○○○○ ○○ ○○○○○ ○○ ○○○ ○○○ ○○○○ ○○○ ○ ○ ○○○○○ ○○○ ○○ ○○○○ ○○ ○○ ○ ○○○○ ○○○○○ ○ ○○○ ○○○○○○ ○○ ○○○○ ○○ ○○○ ○ ○○○ ○○○ ○○○○○ ○○○○○○ ○○○○ ���� ���� ���� �������������������� ���(���� ���� ����� )���●���� ○����Figure 8: Comparison of mean cell-speci c coecients for di erent estimation precedures. (a) ASK14 base model; (b) BSSA14 base model; (c) CB14 base model. (�) ��� ��� ��� ��������������������� ���������� ��������(°)����������� (�) ��� ��� ��� ����������������������� ���������� ��������(°)�������������� �� ��������� ������ �� ��������� �%/��%�������� �� ��������� ���� ��������� ������ ��������� Figure 9: (a) Correlation of spatially varying event terms, depending on separation distance. Shown are the range of correlations basd on the posterior distribution of `eq, as well as correla- tions based on point estimates of the posterior (mean/median of posterior distribution of `eq). (b) Covariance of spatiall varying event terms, taking into account uncertainty of eq. The ergodic base model is ASK, the nonergodic adjustment model is NNGP. models, we used a truncated normal distribution (cf. Equation (12)) as the prior distribution for the cell-speci c e ects, while we use an untruncated normal distribution in INLA. The number of negative cell-speci c attenuation coecients is small, and in a forward prediction they should be set to zero. The posterior distributions of the hyperparameters `eqand eq, governing the distribution of the spatially varying even terms, are shown in Figure 5. The impact of the uncertainty in the hyperparameters on the correlation/covariance for di erent separation distances j~te~t0 ej is shown in Figure 9, for the ASK14 NNGP model. For each sample of the posterior distri- bution, we calculate the correlation and covariance for di erent separation distances according to Equation (8). The 5%,95% quantiles, as well as mean and median, of the resulting correla- tions/covariances are shown in Figure 9. We also show the correlations/covariances, calculated with the mean and median of the posterior distributions. The range of possible correlation/covariance values is quite broad, re ecting the uncertainty in the values of `eqand eq. A full nonergodic PSHA would take the uncertainty in the correlation/covariance into account. However, we will show in Section Prediction of Spatial Terms that the impact of his uncertainty on prediction of spatial terms is very small. 14 ���� ��������� �����������������(�)����� ������ ����Figure 10: Runtime in hours for the di erent models. For the Stan models, the average time of four chains is plotted. 4.1 Runtime of Di erent Models Figure 10 shows the absolute runtime of the di erent models in hours. We can see that using the NNGP approximation in Stan results in a ve-fold speedup compared to the full Stan model. Since the estimated parameters between these two models are almost indistinguishable, this means that the use of the NNGP model can be used to estimate a Bayesian nonergodic VCM for larger data sets and/or more target variables. However, the runtime for the NNGP Stan model is sill quite long with about 20h. By contrast, the INLA model takes only about 1-2h to run, with very similar results to the other models. One should not overinterpret absolute runtimes, as they depend on the load of the computer. We also ran 1000 warm-up and post warm-up samples for the Stan models; in general, the adaptation phase of HMC is quite fast, so one could potentially reduce the number of the warm-up iterations, which would lead to a smaller runtime. 5 Calculating Predictions Following Landwehr et al. (2016), we have cast the nonergodic adjustment models as VCMs with GP priors on the adjustment terms. Hence, the predictive distribution of adjustment termcfor a new location ~tis a normal distribution with predictive mean and predictive standard deviation . we can calculate the predictive distribution by conditioning on observed values (cf. Equations (2.18) and (2.19) in Rasmussen and Williams (2006)) p(cj~c) =N(; ) (32) =~kTK1~c (33) 2 =k~kTK1~k (34) whereKis the covariance matrix of the observed systematic e ects, and ~kis a vector whose components are the covariance of the new location ~twith the observed locations ~t. The entries ofKand~kcan be calculated from the covariance function with the appropriate coordinates. In our case, cis either the event or site adjustment, associated with the coordinates of the hypocenter or the site, respectively. Note that the formulation in Equation (32) is di erent from the corresponding Equation in Landwehr et al. (2016). They integrated out all the nonergodic coecients, and thus condi- 15 tioned on the observed PSA values to calculate the predictive distribution. This means that the matrix Kin Equations (33) and (34) needs to include the aleatory variance, since the ob- servations are not noise-free (cf. Equation (2.20) to (2.21) in (Rasmussen and Williams, 2006)). By contrast, in our case we explicitly estimate the adjustment terms ~cas latent variables; they are not integrated out. Since the latent adjustment variables are theoretically considered noise-free, the matrix Kdoes not include aleatory variability terms. The values of ~c, however, are estimated from data, and thus are associated with epistemic uncertainty. If this uncertainty is ignored, then the estimated predictive uncertainty will potentially be too small. In general, one needs to integrate over the posterior distribution of the parameters governing the nonergodic adjustments to account for their uncertainty in the predictive distribution. Thus, the posterior predictive distribution of adjustment terms can be calculated by p(cjD) =Z p cj~c; ;` p ~c; ;`jD d ~c; ;` (35) Given samples from the posterior distribution, the integral can be approximated bu p(cjD)1 MMX m=1p cj~c(m); (m);`(m) (36) where superscript ( m) denotes the mth sample from the posterior distribution. Calculating the predictive distribution p(cjD) from samples of the posterior distribution requires calculating and inverting the covariance matrix Kfor each sample of the posterior distribution, which can be computationally expensive. If one considers only uncertainty in the latent variables ~c, and assumes that this uncertainty is Gaussian, then one can calculate the predictive distribution as (Lavrentiadis et al., 2021) =~kTK1~ c (37) 2 =k~kTK1~k+~kTK1 c(K1~k)T(38) where~ cis a vector comprising the means of the posterior distributions of the estimated ~c values, and cis a diagonal matrix whose entries are the variances of the posterior distributions of the estimated ~cvalues. For the Stan models, one can calculate the predictive distribution at a new location according to Equations (37) and (38) (using a point estimate of the hyperparameters to calculate the covariances k,~k, andK), or via Equation (36) if one wants to take the full uncertainty into account. The R-INLA package provides tools to sample from the posterior distribution of the model at the mesh nodes, which can then be interpolated to a new location. The samples from the INLA model at new locations implicitly take the full uncertainty of all locations into account. 5.1 Prediction of Spatial Terms Figure 11 shows maps of the mean of systematic event terms ceqover California, for the ASK14 base model. We calculate the values for the mesh nodes shown in Figure 2. For the INLA model, the output contains mean and quantiles of the marginal posterior distribution for the systematic terms at each mesh node. For the Stan models, we calculate the values of the systematic terms at each mesh node based on Equation (37) and (38); for the values of the hyperparameters eqand`eq, we use the median of the respective posterior distributions. As before, Figure 11 shows relatively small di erences between the Stan models and the INLA model. This is concordant with Figure 5, where the posterior distributions of eqand 16 (�) (�) (�) -���-��� � ��� ���Figure 11: Spatially varying event terms ceq, for di erent estimation methods, and ASK14 as the base model. (�) (�) (�) � ���� ���� ���� Figure 12: Epistemic uncertainty associated with spatially varying event terms ceq, for di erent estimation methods, and ASK14 as the base model. `eqare similar. We show similar maps for the other ergodic base models, as well as maps for the systematic site terms cstat, in the electronic supplement. The corresponding epistemic uncertainties are shown in Figure 12. The maps are quite similar, with INLA predicting slightly higher uncertainty. The uncertainty associated with the spatial event adjustment terms ceqranges from 0.075 in the vicinity of observed events to 0.131 in regions where observed events are sparse. In Figures 13 and 14 we compare the predictive mean and standard deviation of the Stan ASK14 NNGP model, calculated for each mesh node, and based on di erent methods to take into account the epistemic uncertainty of the latent function values ~ceqand hyperparam- eters eqand`eq. We compare four di erent variations: •Model A: Taking into account the epistemic uncertainty of the latent function values ~c, and the hyperparameters eqand`eq, based on posterior samples. The predictive distributions at each mesh node are calculated according to Equation (36). •Model B: Similar to Model A, but now we only take into account epistemic uncertainty in the late function values ~ceq. For the hyperparameters, we use the median of the posterior distribution. Predictions are again calculated according to Equation (36). 17 (�) (�) -�����-����� � ����� �����Figure 13: Di erence in mean predictions of spatially varying event terms, for the ASK14 NNGP model, based on di erent methods to take into account uncertainty. (a) mean prediction of model B minus mean prediction of model A. (b) mean prediction of model C minus mean predictionof model A. See text for description of models A, B, C. •Model C: The predictive distribution is calculated according to Equations (37) and (38). For the hyperparameters, we use the median of the posterior distributions. The matrix cis a diagonal matrix whose entries are the variances of the posterior distributions of the estimated ~cvalues. •Model D: Similar to Model B, but the cis the full covariance matrix of the estimated ~c. We use model A as the reference model, since it includes uncertainties on all parameters, and compare the results from the other three models against the predictions from model A. The predictive means for model C and D are the same. The comparisons are shown in Figure 13 for the predictive mean and Figure 14 for the predictive standard deviation , for all mesh nodes. In general, the di erences in the predictive means (Figure 13) are quite small, in particular when compared to the overall range of the nonergdic adjustments (Figure 11). The standard deviations for model A and model B are very similar, which implies that the uncertainty in the hyperparameters eqand`eqhas only a small in uence on the overall epistemic uncertainty of the predictions. There are larger di erences in the standard deviations of model C, which shows lower values of close to the observed data. Overall, these di erences are still small compared to the overall predictive uncertainty. Taking into account the covariance between the estimated latent function values (model D) brings the values of closer to the ones of model A, as seen in Figure 14 (c). It is important to note that the computation of the predictive distribution according to Equation (36) (model A and B) requires repeated calculations of ;kand ;k, for di erent samples kof the posterior distribution. By contrast, model C and D require only one calculation step, and are thus more ecient computationally. 6 Discussion With an increasing number of ground-motion recordings, it is possible to relax the ergodic as- sumption in GMM modeling, and consider systematic, repeatable source, site, and path e ects. This provides challenges, both from a conceptual point of view (how should the systematic e ects be modeled), as well as computational. Landwehr et al. (2016) introduced a VCM to 18 (�) (�) (�) -�����-����� � ����� �����Figure 14: Di erence in standard deviations of predictions of spatially varying event terms, for the ASK14 NNGP model, based on di erent methods to take into account uncertainty. (a) standard deviation of prediction of model B minus of model A. (b) standard deviation of prediction of model C minus of model A. (b) standard deviation of prediction of model D minus of model A. See text for description of models A, B, C, D. model nonergodic GMMs, based on GPs. While conceptually attractive due to their simplicity, as well as the ease to estimate predictive uncertainties, GPs pose diculties when applied to large data sets, since they do not scale well to large data sets. Here, we have investigated di erent methods to t Bayesian nonergodic GMMs, based on INLA (Rue et al., 2009) and MCMC sampling via Stan Carpenter et al. (2017). For Stan, we t both a full (without any approximation) as well as an approximated model based to NNGPs (Datta et al., 2016b). Since the development of full nonergodic GMMs is still in its infancy, it is important to investigate di erent methods to estimate them. Though the INLA and Stan models are based on slightly di erent assumptions, we nd that the results are quite similar between the models. In terms of runtime of the regression models, INLA provides a clear advantage in our ex- periments (cf. Figure 10). We also nd that the NNGP approximation has a clear speed-up over over the full Bayesian MCMC model without any approximation. Further avenues of investigation can include other approximations, such as combination of varational inference and MCMC (Margossian et al., 2020) or reduced-rank GP regression (Solin and S arkk a, 2020; Riutort-Mayol et al., 2020). Recently, Dambon et al. (2021) proposed a method based regu- larized maximum likelihood estimation to estimate VCMs with many varying coecients on large data sets, which can be explored as well. As an alternative to nonergodic GMMs based on GP VCMs, Caramenti et al. (2020) proposed a method based on geographically weighted regression (GWR). Here, we are mainly concerned with the comparison of di erent estimation approaches for nonergodic GMMs. Both methods are Bayesian models. We believe that Bayesian models provide advantages, as they allow to incorporate prior information in a probabilistic way, and can be easily updated when new data becomes available (Sta ord, 2019). We have cast the nonergodic models as adjustments to ergodic base models, but both INLA and Stan allow to easily combine the tting of on-spatially varying and spatially varying coecients. As shown in Kuehn et al. (2019) (and also seen in Figure 7), the average anelastic attenuation coecient changes when estimating cell-speci c path e ects, so taking into account nonergodic e ects could have an e ect on other parameters as well. It is easy to extend the nonergodic models, both for INLA and for Stan, to incorporate other spatially varying terms. For example, Landwehr et al. (2016); Caramenti et al. (2020); 19 Sung et al. (2021) account for a spatially varying VS30-scaling in their models. One could also incorporate di erent spatial length scales for the same parameter, to account for regional and local variation. In this case, one should use informative prior distributions to avoid problems with identi ability of the di erent length scales. In our models, the spatially varying adjustment terms are modeled with a stationary and isotropic correlation function (cf. Equations (8) and (25)). There is evidence that spatial corre- lations of ground-motion within-event are better modeled with a non-stationary or anisotropic correlation function (Kuehn and Abrahamson, 2020; Abbasnejadfard et al., 2020), and one would expect this to be true for spatially varying coecients as well. INLA models have been developed to account for spatial non-stationarity (Ingebrigtsen et al., 2014, 2015; Bakka et al., 2019); it should be investigated whether incorporating non-stationarity into the nonergodic models leads to improved performance. One important aspect when relaxing the ergodic assumption is to take into account epistemic uncertainty. The nonergodic adjustments are by nature estimated from smaller amounts of data, and are thus less well constrained compared to ergodic GMMs. In a VCM, there is parameter uncertainty (cf. Figure 9), as well as uncertainty for predictions at new locations. Figures 13 and 14 indicate that uncertainty in the hyperparameters only has a small in uence on the overall predictive uncertainty. Acknowledgements We used the python package ArviZ (Kumar et al., 2019), as well as R-package posterior (B urkner et al., 2014) to perform analysis of the Stan models. This work was sponsored by the PG&E Geosciences Department. We would like to thank the nonergodic working group at UC Berke- ley/UCLA for discussions, in particular Jennifer Donahue and Kathryn Woodell.