A Bayesian latent factor to estimate multivariate empirical ground-motion models is proposed. In the factor model, event terms, site terms, and within-event/within-site residuals for each response oscillator period are modeled as linear combinations of a small number of latent factors, thus introducing inter-period correlation in the model. The factor model provides a lower-dimensional approximation to the full residual correlation structure. The factor model is demonstrated on Italian ground-motion data. Results indicate that the factor model leads to similar behavior as a traditional ground-motion model where each target is estimated separately. 1 Introduction Empirical ground-motion models (GMMs) are an indispensable tool for probabilistic seismic hazard analysis (PSHA). A GMM provides an estimate of the distribution of a ground-motion parameter (such as peak ground acceleration (PGA) or the response spectrum (PSA)), con- ditioned on a particular source/site scenario. Empirical GMMs are typically estimated via regression, which is usually carried out separately for each target variable. Correlations be- tween response specral periods are useful for ge.g. generating conditional mean spectra ( Baker , 2011) or joint hazard ( Bazzurro and Cornell , 2002). Correlation models are typically estimated from residuals (e.g. Baker and Bradley , 2017; Baker and Jayaram , 2008; Bradley , 2011; Huang and Galasso , 2019; Kotha et al. , 2017) to empirical GMMs. Some studies (e.g. Arroyo and Ordaz , 2010a,b; Kuehn and Scherbaum , 2015) have proposed methods to estimate multivariate GMMs, where all ground-motion parameters are regressed simultaneously together with their full covariance. Modeling the full spectrum during regression can be useful, such as when taking into account uncertainty in predictor variables ( Kuehn and Abrahamson , 2018) or for truncated regression ( Kuehn et al. , 2020b) due to triggering. Hee, we model the residual structure in empirical GMMs with a factor model, which accounts for inter-period correlation by modeling the relation between periods via a small number of ∗kuehn@ucla.edu 1 latent factors. Factor models are widely used to model the covariance structure in multivariate data via dimensionality reduction, dating back to Spearman (1904) and Thurstone (1935). In the next section, we describe the general model and how it can be use to model the covariance structure of event terms, site terms, and residuals in empirical GMMs. We then estimate the factor model GMM on Italian ground-motion data, and compare the results with a traditional GMM which is estimated separately for each target variable. 2 Factor Model GMM The goal of an empirical GMM is to model the distribution of several ground-motion parameters of interest, given earthquake source and site related parameters such as magnitude, distance, or time-averaged shear wave velocity in the upper 30m, VS30. In the following, let us assume a number of Pground-motion target variables. Typically, empirical GMMs are estimated separately for each target variable. Each target is usually modeled as lognormally distributed, which means that its logarithm is normally distributed: µp=fGMM (⃗ cp;⃗ x) +δBp+δSp yp∼N(µp, ϕSS,p)(1) where pis an index for the target variable, and ypis the logarithm of the pth target variable. ⃗ x is a vector of predictor variables (e.g. magnitude, distance, VS30),⃗ cpis a vector of coefficients for each target variable, and fGMM (⃗ cp;⃗ x) describes the functional dependence of the target variable on the predictor variables. δBandδSare event and site term, respectively, and describe a constant offset for each record from individual events/stations ( Al-Atik et al. , 2010). ϕSSis the standard deviation of the within-even/within-site aleatory variability. The event terms and site terms are typically modeled as independent random variables, normally distributed with mean zero and standard deviations τpandϕS2S,p, respectively ( Al-Atik et al. , 2010). One can write the model as as a multivariate normal distribution ⃗ y∼MN(⃗ µ,ΦSS) (2) where Φ SSis the covariance matrix between the different target variables. When estimating empirical GMMs, the models are typically estimated separately for each target variable, which means that Φ SSis a diagonal matrix with entries Φ SS,pp =ϕ2 SS,p. Similarly, event and site terms across periods can be modeled as multivariate normal distributions with a mean vector of zeros and diagonal covariance matrices Tand Φ S2S. One can relax the assumption of independence, and estimate the joint distribution of all target variables during model fitting. This requires to perform a regression for all target periods at the same time, with estimation of full covariance matrices Φ SS,T, and Φ S2S. The truncated regression model proposed in Kuehn et al. (2020b) relies on modeling the joint distribution of different target variables. Previous studies of multivariate GMM regressions (e.g. Arroyo and Ordaz , 2010a,b; Kuehn and Abrahamson , 2018; Kuehn et al. , 2020b) modeled the covariance matrices as full-rank, which can be computationally expensive. Here, we propose a GMM factor model, where structure across periods is introduced by a linear combination of latent (unobserved) factors. The factors are of lower dimension than the targets. 2 In the following, we briefly describe the factor model, following the exposition given in Papastamoulis and Ntzoufras (2022). In a factor model, the target variables are modeled as a linear combination of the latent factors ⃗ y=⃗ µ+ Λ⃗ η+⃗ ϵ (3) where Λ is a P×Dfactor loadings matrix, ⃗ ηareDlatent factors, an ⃗ ϵarePindependent errors, typically called idiosyncratic errors . For D < P , this is an approximation of the full problem with a lower-dimensional problem. The errors are independent from the factors, and independent normal random variables ⃗ ϵ∼N(0,Σ) (4) where Σ is a diagonal matrix. The distribution of ⃗ yconditional on the latent factors is ⃗ y|⃗ η∼N(⃗ µ+ Λ⃗ η,Σ) (5) and the marginal distribution is ⃗ y∼N(⃗ µ,ΛΛT+ Σ) (6) The term ΛΛTdescribes the shared variation between targets ( Heaps , 2022). Each factor is assumed to be independent and normally distributed with standard deviation one λp∼N(0,1) (7) For the factor GMM, we model event terms, station terms, and within-event/within-station residuals across multiple periods with a factor model. For event terms and station terms, the model is⃗δB= Λ E⃗ ηE ⃗δS= Λ S⃗ ηS(8) ηEandηSare the latent factors for events and stations, and Λ Eand Λ Sare the associated factor loadings matrices. We omit the idiosyncratic errors for event and site terms, since we found that their values were very small, and they hindered model estimation. The standard deviations of the event and site terms for each period can be calculated as ⃗ τ=q diag(ΛEΛT E) ⃗ϕS2S=q diag(ΛSΛT S)(9) i.e. as the square-root of the diagonal of the shared variation matrix. The within-event/within-site residuals are modeled as ⃗δWS = Λ R⃗ ηR+⃗ ϵR (10) where subscript Rstands for record, and the other terms are defined as before. The idiosyncratic errors ⃗ ϵare each normally distributed with mean zero and standard deviation ϕ0,p. Thus the total within-event standard deviation can be calculated as ⃗ϕSS=q ⃗ϕ2 λ+⃗ϕ2 0 ⃗ϕ2 λ=diag(ΛRΛT R)(11) 3 2.1 Implementation The parameters of the models are estimated with Bayesian inference ( Gelman et al. , 2013) via Markov Chain Monte Carlo (MCMC) sampling. The models are implemented in the proba- bilistic programming language Stan ( Carpenter et al. , 2017; Stan Development Team , 2022) and run from the computing environment R(R Core Team , 2022) using the package cmdstanr (Gabry and ˇCeˇ snovar , 2021). Bayesian inference of factor models with MCMC is not straightforward ( Papastamoulis and Ntzoufras , 2022). The factor model is not uniquely identified, as any rotation of the factor loadings matrix and the factors leads to the same model. To avoid problems with rotation, the factor loadings matrix is often constrained such that λij= 0 for j > i andλii>0 (Geweke and Zhou , 1996; Lopes , 2014; Lopes and West , 2004); thus, the upper triangle is zero. Here, we order the target variables from short periods to long periods, and constrain the factor loadings matrix such that the lower triangle s zero. The reason is that at longer periods, there is less data due to filtering, which makes the association to the factor loadings harder to detect. We illustrate the model using a slightly modified version of the model and data of the ITA18 GMM Lanzano et al. (2019a). The model is simple (linear in coefficients), and the number of data points is not too high. We illustrate the model using 18 periods (mainly a computational reason) and eight factors. For comparison, we also fit a Bayesian model where each target variable is estimated separately. We run two chains with 500 warm-up samples and 500 estimation samples each, resulting in 1000 posterior samples. 3 ITA18 Model and Data The ITA18 GMM ( Lanzano et al. , 2019a) is a model developed for Italian ground-motion data. The data is taken from the Engineering Strong Motion (ESM) flatfile Lanzano et al. (2019b), enhanced by some small regional Italian records from the ITalian ACcelerometric Archive (ITACA) database ( Luzi et al. , 2008, 2017; Pacor et al. , 2011), as well as 12 worldwide events in the magnitude range 6.1 to 8. In total, ITA18 is based on 5607 records from 146 events and 1657 stations. Here, we use the subset of ITA18 from the ESM data base, which consists of 4541 Italian records from 137 events and 886 stations. The magnitude-distance scatterplot is shown in Figure 1. The functional form of ITA18 relating the target variable to the predictors is y=a+b1(MW−Mh) 1(Mw≤Mh)+b2(MW−Mh) 1(Mw>Mh) + [c2+c1(MW−Mh)] log10q R2 JB+h2+c3q R2 JB+h2 +k log10VS30 8001(VS30≤1500)+1500 8001(VS30>1500) +δB+δS+δWS(12) where yis the target variable (in log10 units), and a,b1,b2,c1,c2,c3, and kare coefficients 4 4567 0.1 1.0 10.0 100.0 RJB (km)MFigure 1: Magnitude-distance scatterplot of the data. 170000175000180000185000190000 2.5 5.0 7.5 10.0 12.5 Dlp Figure 2: MAP log density versus number of factors, estimated for within-event/within-station residuals. that are estimated during the regression. The predictor variables are moment magnitude MW, Joyner-Boore distance RJB,VS30. The original ITA18 model includes two coefficients describing the dependence of yon fault type; we omitted these parameters Mh,Mref, and hare parameters that are fixed before model estimation. ITA18 was developed for PGA, PGV, and 36 response spectral periods from T= 0.01s toT= 10 .s. In this work, we estimate the different models for PGA and periods from 17 periods from T= 0.01−3, i.e. P= 18. A longer periods, data drops out due to filtering. We use D= 8 for the umber of latent factors. The number is chosen as follows: we take the within-event/within-station residuals of the separate model, and estimate a simple factor model like Equation (3) via maximum-a-posteriori estimation for different numbers of factors. Figure 2 shows the value of the maximised log-density (lp) versus the number of factors. At about D= 8, the log-density reaches a plateau. For the factor model, we also provide reasonable starting values for the MCMC algorithm, as this can to avoid label switching in the entries of the factor loadings matrix ( Deschamps , 2018). For the coefficients, the starting values are taken from the posterior distribution of the separate model, while stating values for the the factor loadings are generated by randomizing the values obtained from optimization based on within-event/within-station residuals. Starting values for the latent factors are generated by sampling from the prior (a standard normal distribution). 5 0.00.10.20.3 0.01 0.10 1.00 10.00 T (s)φS2SFactor Separate 0.00.10.20.3 0.01 0.10 1.00 10.00 T (s)τ 0.00.10.20.3 0.01 0.10 1.00 10.00 T (s)τ 0.200.250.300.350.40 0.01 0.10 1.00 T (s)σσSSσTFigure 3: Estimated standard deviations vs period. Plotted is the mean of the posterior distribution, error bars indicate 5% to 95% intervals of the posterior. The periods have been slightly moved to avoid plotting points at the same place. The blue line shows the value of ITA18. 4 Comparison between Factor and Separate Model In his section, we show comparisons of the estimated factor model GMM as well as the model where each target variable is estimated separately. For the factor model, there are some warn- ings that the tree-depth of the NUTS algorithm ( Hoffman and Gelman , 2014) has reached the maximum limit, and some values of the ˆR-statistic ( Vehtari et al. , 2021) are larger than 1.05 for some parameters. Thus, the results for the factor model should be interpreted with care, but we believe that the estimated values are reasonable. Figure 3 shows the estimated standard deviations ϕS2S,τ,ϕ. In addition, we show total standard deviation σ2 T=ϕ2 S2S+τ2+ϕ2 SS, and single-station standard deviation σ2 SS=τ2+ϕ2 SS. In general, these are very similar between the models. The factor model leads to some slightly larger values. Due to the modeled correlation, the individual event/site terms and residuals have to obey some additional constraints. Figure 4 shows the values of the components of the within-event standard deviation (cf. Equation (11)) for the full factor model. Most of the within-event/within-station variability is explained by the shared variation, with a larger share of idiosyncratic errors at longer periods. The estimated coefficients are shown in Figure 5. In general, the coefficients are very 6 0.000.050.100.150.200.25 0.01 0.10 1.00 T (S)φSSφ φλ φ0 0.000.050.100.150.200.25 0.01 0.10 1.00 T (S)φSSφ φλ φ0Figure 4: ϕSSand its components for the factor model. similar between models. The coefficient b2, which controls the large magnitude scaling, shows the largest difference. The reason for this is probably that there are not may large events, and the modeled correlation between event terms enforces some consistent behavior across periods. Figure 6 shows the empirical correlation between the posterior samples of some periods across periods, i.e. we calculate the correlation at two periods T1andT2 ρ(aT1, aT2) (13) in this case for the intercept a. Figure 6 shows the correlations with respect to T= 0.5s, for the intercept a, the large magnitude scaling b2, and the anelastic attenuation c3. While the priors for each coefficient are independent for both the separate and factor model, there is strong correlation between the coefficients at different periods in the factor model. By contrast, the correlation is (expected) close to zero for the separate model. This would have a consequence if one calculates joint predictions for several periods and wants to take the epistemic uncertainty in the correlations into account. Figure 7 shows correlations between different periods, calculated from the shared variation ΛΛT, for event terms, site terms, and within-event/within-site residuals. For comparison, we also calculate empirical correlations from the estimated terms for the separate model, and show the correlation model of Baker and Jayaram (2008) (note that this model is for correlation between total residuals). The correlations are all quite similar, which is in line with previous research that has found that period-o-period correlations of the response spectrum are stable across different regions and ranges of predictor variables ( Baker and Bradley , 2017). 5 Discussion and Conclusions We have presented a factor model GMM, to estimate empirical GMMs for multiple target variables at the same time, while also accounting for inter-period correlation in event terms, site terms, and within-event/within-site residuals. The factor model GMM is a low-rank ap- proximation to a full covariance matrix, and thus models the residual structure using a lower- dimensional space. In general, results from the factor model GMM and a separate model (i.e. all targets are estimated without accounting for correlation) are similar, both in terms of the estimated coeffi- 7 234 0.01 0.10 1.00 10.00 T (s)a Factor Separate 0.000.250.500.75 0.01 0.10 1.00 10.00 T (s)b1 −0.250.000.250.500.75 0.01 0.10 1.00 10.00 T (s)b2 0.10.20.3 0.01 0.10 1.00 10.00 T (s)c1 −1.5−1.4−1.3−1.2 0.01 0.10 1.00 10.00 T (s)c2 −0.004−0.0020.000 0.01 0.10 1.00 10.00 T (s)c3 −0.8−0.6−0.4 0.01 0.10 1.00 10.00 T (s)kFigure 5: Estimated coefficients vs. period. Dots show the mean of the posterior distribution, the uncertainty interval comprises the 5% and 95% quantiles of the posterior distribution. The blue line shows the coefficients of ITA18. 0.000.250.500.751.00 0.01 0.10 1.00 T (s)ρ(a) 0.000.250.500.751.00 0.01 0.10 1.00 T (s)ρ(b2) 0.000.250.500.751.00 0.01 0.10 1.00 T (s)ρ(c3) Figure 6: Empirical correlation between constant ( a), large magnitude scaling ( b2), and anelastic attenuation ( c3) coefficients, conditioned on T= 0.2s, estimated from the posterior samples. 8 0.000.250.500.751.00 0.01 0.10 1.00 T (s)ρE Factor full Separate BJ08PGA 0.000.250.500.751.00 0.01 0.10 1.00 T (s)ρET = 0.200s 0.000.250.500.751.00 0.01 0.10 1.00 T (s)ρET = 1.000s 0.000.250.500.751.00 0.01 0.10 1.00 T (s)ρSPGA 0.000.250.500.751.00 0.01 0.10 1.00 T (s)ρST = 0.200s 0.000.250.500.751.00 0.01 0.10 1.00 T (s)ρST = 1.000s 0.000.250.500.751.00 0.01 0.10 1.00 T (s)ρRPGA 0.000.250.500.751.00 0.01 0.10 1.00 T (s)ρRT = 0.200s 0.000.250.500.751.00 0.01 0.10 1.00 T (s)ρRT = 1.000sFigure 7: Estimated correlations, and comparisons with model of Baker and Jayaram (2008). Dashed indicates correlations are estimated from residuals. 9 cients and standard deviations. This raises the question whether it is useful to account for the inter-period correlation structure during regression. While results are quite similar, there are use cases, such as truncated regression ( Kuehn et al. , 2020b), or accounting for uncertainty in predictor variables ( Kuehn and Abrahamson , 2018), which require modeling the full spectrum. In the factor model GMM, we account for correlations in the residuals, while the coefficients are modeled separately. The factor model induces correlation in the posterior distribution of the coefficients (cf, Figure 6), but no structure or smoothing is applied for the behavior of the periods across periods. One way forward could be to impose a smooth relation of the coefficients. Bortolotti et al. (2022) describe the use of functional data analysis to obtain smooth coefficients. One could also model the functional dependency of the coefficients on periods using a Gaussian process, as in Kuehn et al. (2020a). The factor model requires choosing the number of latent factors. Often, the maximm number of factors is taen t be the Lederman bound ϕ(P) =2P+ 1−p (8P+ 1) s(14) since Σ (the covariance matrix of the idiosyncratic errors) is globally identifiable when P < ϕ(P) (Bekker and Ten Berge , 1997; Heaps , 2022; Papastamoulis and Ntzoufras , 2022). While there are methods proposed to identify the number of latent factors during model estimation (Fruehwirth-Schnatter and Lopes , 2018; Heaps , 2022), here we fixed the number based on the maximum estimted log density when optimizing for a smaller problem (cf. Figure 2) for simplicity. The components of variability and uncertainty in GMMs ( τ,ϕS2S,ϕSS) are sometimes modeled as dependent n predictor variables; in particular he between-event variability τis often assumed to decrease for large magnitudes. In the factor model, the value of τis the square root of the diagonal of Λ EΛT E. Thus, dependence of τon magnitude needs to be modeled through dependence on the factor loadings matrix, which would lead to a magnitude-dependent correlation structure. Alternatively, one could model the distribution of the latent factors as magnitude dependent. Here, we have used the factor model to estimate an empirical GMM, but there are other applications in the context of ground-motion modeling and hazard/risk analysis. The latent factors for events and sites ( ηEandηS) are assumd to be independent, but they could be modeled as spatially correlated ( Qiu et al. , 2018; Thorson et al. , 2015a,b). This opens an avenue for modeling multi-variate nonergodic GMMs ( Lavrentiadis et al. , 2022). Similarly, a spatial factor model could be used for mdelig cross-correlations of ground-motion parametes, which is useful for regional seismic risk analysis (.g. Miller and Baker , 2015). Such a model provides a lower dimensional approximation to the full cross-covariance, which is similar in spirit to the model of Markhvida et al. (2018). In general, we believe that the factor model provides a convenient tool for modeling the full response spectrum in empirical GMMs in case this is needed. 6 Data and Code Data is available from ( Lanzano et al. , 2019b) and Lanzano et al. (2019a). Stan and R code to run the factor and separate models is available at https://github.com/nikuehn/ 10 FactorModelGMM . Acknowledgements The regressions are carried out in R(R Core Team , 2022) usig package cmdstanr (Gabry and ˇCeˇ snovar , 2021). Plots are made with ggplot2 (Wickham , 2016), and data organization is done using tidyverse (Wickham et al. , 2019) and posterior (B¨ urkner et al. , 2022). Partial support of Pacific Gas & Electric Company and California Department of Trans- portation are gratefully acknowledged. Any opinions, findings, and conclusions or recommen- dations expressed in this material are those of the author and do not necessarily reflect those of the sponsors. 11