We provide a simple introduction to the estimation of ground-motion models via Bayesian inference and the probabilistic programming language Stan. We show one ca implement a simple ground-motion model in Stan, and how can run the program from the computer environment R. We also show how one can access the results, and plot summaries of estimated parameters. A large number of di erent Stan models for the development of GMMs is available at https://github.com/pstafford/StanGMMTutorial . 1 Introduction Empirical ground-motion models (GMMs) are an essential ingredient for probabilistic seismic hazard analysis (PSHA). They are typically estimated from data via regression, often using maximum-likelihood. Recently, Sta ord (2019) provided an argument that GMMs should be estimated via Bayesian inference, as it provides a toolchain for integrating new data into ex- isting GMMs. Bayesian methods have been used for estimation of a GMM a couple of times (e.g. Arroyo and Ordaz, 2010,a; Kowsari et al., 2019, 2020; Kuehn and Abrahamson, 2018; Kuehn et al., 2019; Kuehn and Scherbaum, 2015, 2016; Ordaz et al., 1994; Rahpeyma et al., ∗kuehn@ucla.edu †p.sta ord@imperial.ac.uk 1 2018; Sta ord, 2014, 2019; Wang and Takada, 2009), but have not found widespread adoption, which might be due to the fact that Bayesian methods are sometimes viewed as complicated. Here, we want to give a short tutorial how to implement Bayesian regression models in Stan (Carpenter et al., 2017) for the estimation of GMPEs. Stan is a program that implements Markov Chain Monte Carlo inference to estimate the posterior distribution (similar programs are WinBUGS/OpenBUGS (Lunn et al., 2009, 2000) or JAGS (Plummer, 2003), though these are using Gibbs sampling, as opposed to Hamiltonian Monte Carlo (HMC) like Stan). We do not want to go into philosophical arguments about Bayesian statistics; we treat Bayesian infer- ence as a practical tool that helps to combine prior knowledge and complex models. However, we note that GMM development typically requires a lot of judgement and prior information, so that the model behaves well and extrapolates well to large magnitudes and sparse distances (data ranges which are typically sparse). Here, prior distributions are helpful to constrain the model coecients to reasonable values. Our goal in this paper is very speci c: to provide a collection of models, coded in Stan, that are useful for various aspects of GMM development. This is not a tutorial or introduction to Bayesian statistics. It is also not a tutorial to MCMC. We provide some information about these topics; in particular, model assessment is important in the context of MCMC. However, we recommend that the interested reader delves deeper into the topic. There are textbooks on Bayesian statistics that include Stan code (e.g. Gelman et al., 2013; Kruschke, 2015; McElreath, 2020). For thoughts on a principled Bayesian work ow, see Gelman et al. (2020). As stated before, in this work we focus on the implementation of GMMs in Stan. We provide a large collection of models at https://github.com/pstafford/StanGMMTutorial . Example Models include truncated regression models (Kuehn and Abrahamson, 2020; Kuehn et al., 2020b), models with correlated target variables (Kuehn and Scherbaum, 2015), measurement uncertainty models (Kuehn and Abrahamson, 2018; Sta ord, 2014), partially nonergodic models (Kuehn et al., 2020a; Kuehn and Scherbaum, 2016; Sta ord, 2014), and models including spatial correlation Huang and Galasso (2019); Jayaram and Baker (2009); Kuehn and Abrahamson (2020); Sta ord et al. (2019). We also show one can code models with heteroscedastic variances, or implement Bayesian robust regression (Kuehn et al., 2020a). 2 2 Bayesian Inference Broadly speaking, the goal of statistical inference is to estimate the values of parameters of a statistical model from data. Given a data set D=h X;~Yi , whereXis a matrix of predictor variables and ~Yis a vector of target variables, we are interested to infer the parameters of a model of the form y=f(~ x;~ c) +" (1) wheref() is a function that relates the target variable to the predictor variables, ~ care coef- cients of the function, and "is an error term, typically assumed to be normally distributed with mean zero and standard deviation . We will call all parameters of the model ~; in this simple model, the parameters comprise the coecients ~ cand the standard deviation . An important quantity in this regard is the likelihood functionL(~), which is a function of the parameters ~of the model. The value of the likelihood function for a given observation iis de ned as the value of the probability density of the observation, given parameters ~: L(~j(yi;~ xi)) =p(Y=yij~ xi;~) (2) For a set of observations (a data set D=f(y1;~ x1);:::; (yN;~ xN)gT), one can calculate the likelihood of the full data set L(~jD) =p(~Y=~ yj~X;~) (3) where~Xis a matrix containing the predictor variables for all observations. If the observations are independent, the likelihood factorizes, and one can calculate the joint likelihood as L(~jD) =NY i=1L(~j(yi;~ xi)) (4) In a maximum-likelihood (ML) setting, the goal is to nd the set of parameters ~that maximizes the joint likelihood, given the data. Note that, given the de nition of the likelihood function in Equation (2), one actually maximizes the probability of seeing the data, given the parameters. By contrast, in Bayesian inference the goal is to estimate the so-called posterior distribution of the parameters, which is the conditional distribution of the parameters, given 3 the observed data. The posterior distribution can be calculated by Bayes' theorem as p(~jD)/p(Dj~)p(~) (5) wherep(Dj~) is the likelihood, and p(~) is called the prior distribution. The prior distribution p(~) allows one to incorporate prior information about some parameters, or to constrain the behavior of the model; this is important in GMM development, as data is generally sparse at large magnitudes and short distances, so the data is not very informative. The prior distribu- tion can also serve as some sort of regularization to avoid over tting to the data. A simple introduction to Bayesian inference can be found in e.g. Spiegelhalter and Rice (2009). There are also many dedicated textbooks, like Kruschke (2015), Gelman et al. (2013) or McElreath (2020). The posterior distribution p(~jD) contains all relevant information about the parameters, conditioned on the data and taking prior information into account. It is thus a good tool to assess the (epistemic) uncertainties of a model. One can calculate summary statistics such as mean/median or standard deviations of individual parameters, as well as correlations between parameters. With the posterior distribution one can also calculate the posterior predictive distribution1, which encapsulates the full uncertainty in a model prediction. It can be calculated by integrating out the uncertainty of the model parameters. Except for simple models, the posterior distribution p(~jD) cannot be calculated analyt- ically, and one has to resort to approximations. One of the most popular tools to perform Bayesian inference is Markov Chain Monte Carlo (MCMC) (Neal, 1993) sampling, in which one constructs a Markov chain whose stationary distribution is the posterior distribution p(~jD). Other approximation methods are variational inference (VI) (Blei et al., 2017), or integrated nested Laplace approximation (INLA) (Rue et al., 2009). In the following, we focus on the pro- gram Stan (Carpenter et al., 2017, https://mc-stan.org/ ). Stan implements both MCMC and variational inference, but we will restrict ourselves to the sampling based methods. It is also possible to perform maximum-a-posteriori (MAP) inference, i.e. nd the set of parameters that maximizes the log-posterior density. There are a few reasons why we focus on Stan. We do not want to discourage people 1https://en.wikipedia.org/wiki/Posterior_predictive_distribution 4 from using other methods; this work is an encouragement to try out Bayesian methods. We use Stan because it allows us to t models of (almost) arbitrary complexity. Some of these more complex models are described at https://github.com/pstafford/StanGMMTutorial . In addition, since Stan is based on HMC (Betancourt, 2017; Betancourt and Girolami, 2015; Neal, 2011), there are multiple tools to assess the model t that can be used to diagnose problems with the model and/or model t (Betancourt, 2016; Vehtari et al., 2020b). Another advantage of using a program such as Stan is that it is actively developed and has a broad user base. While it is possible to implement their own MCMC scheme (e.g. Rahpeyma et al., 2018), the bene ts of using an established software is that one can focus on the modeling part, not the implementation; in addition, maintenance is less of an issue. A dedicated implementation might be faster for a speci c model, but Stan allows one great exibility to try out di erent models. 3 Bayesian Inference for Ground-Motion Models In this section, we motivate the use of Bayesian inference and Stan for the development of GMMs. A typical GMM has the form of a mixed e ects model, to account for repeated measurements of the same event, or at the same station. This gives rise to a model like the following: Yes=f(~ xes;~ c) +Be+Ss+WS es (6) wheref(~ x;~ c) is a function with coecients ~ c, which also depends on predictor variables ~ x.e andsare indices for the event and station. BeandSsare the event and station term for eventeand station s, respectively, and WS is the remaining residual (tpyically called within- event/within-station residual). The function f(~ xes;~ c) is a model for the median prediction for a generic scenario quanti ed by ~ xes. Due to the presence of event and station terms in Equation (6), the observations ~Yesare not independent, and the joint likelihood is a multivariate normal distribution with covariance , with non-zero o -diagonal elements  ijif recordiandjare from the same event or station. In an ML model, the random e ects BandSare considered errors in the constant coecient, and are integrated out of the likelihood; thus, they are not estimated, but rather predicted 5 given the estimated values for the coecients and  (see e.g. Equation (10) of Abrahamson and Youngs (1992)). By contrast, in a Bayesian model BandSare considered parameters of the model that are estimated. We believe that this is a more natural position for GMM development, since one is interested in the di erence of ground-motions at a particular site from the average prediction. Hence, the e ects BandSconstitute physical parameters about which inference is of interest in itself. For a known event term Band station term S, the observations become independent, and the likelihood factorizes: p(~Y=~ yj~X;~) =NY i=1p(yij~ xi;~ c;B eq(i);Sstat (i); SS) (7) where~Xis a matrix of predictor variables, whose rows consist of the vectors of predictor variables for each record. On the right hand side of Equation (7), we have partitioned the full set of parameters ~of the model into the components; ~consists of the event and station terms, the coecients of the functional form ~ c, and the standard deviation of the within-event/within- station residuals SS(the \noise" term). The subscripts eqandstat are indices that connect theith record to its event and station, respectively. As stated in Equation (5), the posterior distribution is proportional to the product of prior distribution and likelihood. For event terms and station terms, the prior distributions are normal distributions with mean zero and standard deviation and S2S, respectively. Since these terms are typically assumed to be independent, the prior distribution for all event terms can be written as p(~B) =QNE i=1p(Bij0; ) (similarly for station terms S). Thus, we can write the log posterior distribution as lnp(~jD) =NX i=1lnp(yijBeq(i);Sstat (i);~ c; SS) +NEX i=1lnp(Bij ) +NSX i=1lnp(Sj S2S) + lnp(~ c) + lnp( ) + lnp( S2S) + lnp( SS) (8) where the rst line is the likelihood, the second line comprises the prior distributions for event and station terms, and the prior distributions for coecients and standard deviations are on 6 the third line. As Equation (8) shows, calculating the value of the log posterior consists of summing di er- ent log probability values. In Stan, this is how one would specify the model. In a Stan program, the statement target += normal_lpdf(y | 0, 10); increments the log density with the value of the expression on the right hand side; in this case, the logarithm of the value of the prob- ability density function (pdf) of a normal distribution with mean zero and standard deviation 10, for a variable y. The construction of a Stan program consists of several statements of the formtarget += such that they add up to the log density of the full posterior distribution. The function normal_lpdf(...|..., ...) is a built-in function for the logarithm of the PDF of a normal distribution ( _lpdf stands for logarithm of pdf). Most commonly used dis- tributions are implemented in Stan in the same way. One can also increment the log density with the statement y ~ normal(0, 10); . This statement is equivalent to target += normal_lpdf(y | 0, 10); except that the former leaves out constant additive terms. For inference, it is only necessary to calculate the value of the log density up to an additive constant. 4 Stan Implementation of Simple GMM In this section, we provide code for a Stan program of a simple GMM and discuss the imple- mentation. We will use a model with the following functional form as an example f(~ x;~ c) =c1+c2M+c3(8M)2+c4ln(R+h) +c5Mln(R+h) +c6R+c7lnVS30 400(9) If the nite-fault term his xed, then the model becomes linear, and can be written as f(~ x;~ c) =~ c~ x (10) with~ x= 1;M;(8M)2;ln(R+h);Mln(R+h);R;lnVS30 400 . For the whole data set, the pre- dictors can be placed into a matrix Xof sizeNK, whereKis the number of (linear) predictors (i.e. K= 7 in this example). The Stan code for the model of Equation (9) is shown below: data { 7 int<lower=1> N; // number of records int<lower=1> N_eq; // number of earthquakes int<lower=1> N_stat; // number of stations int<lower=1> K; // number of predictors matrix[N, K] X; // matrix of predictors vector[N] Y; // target variables int<lower=1, upper=N_eq> idx_eq[N]; // event index for each record int<lower=1, upper=N_stat> idx_stat[N]; // station index for each record } parameters { vector[K] c; // coefficients real<lower=0> phi_SS; // standard deviation for within-event residuals real<lower=0> phi_S2S; // standard deviation of between-event residuals real<lower=0> tau; // standard deviation of site-to-site residuals vector[N_eq] deltaB; // event terms vector[N_stat] deltaS; // station terms } model { // prior distributions c ~ normal(0,10); phi_SS ~ normal(0,0.5); tau ~ normal(0,0.5); phi_S2S ~ normal(0,0.5); 8 deltaB ~ normal(0,tau); deltaS ~ normal(0,phi_S2S); // likelihood Y ~ normal(c * X + deltaB[idx_eq] + deltaS[idx_stat], phi_SS); } The basic structure of the program consists of three blocks, the data {...} block, the parameters {...} block, and the model {...} block. In the data {...} block, the input parameters are de ned. In the case of the GMM, these are the number of data, evens, stations, and linear predictors, the predictor matrix X, the target variable Y, and indices linking records to events and stations. Variables in Stan are typed, meaning that there is a di erence between a real and an integer. An integer variable is de ned as int N; , while a real is de ned as real a; . The target variable is read in as vector[N] Y; , where notation []denotes the size of the vector (similar for matrices). Indexing a vector starts with 1 (i.e., the rst element in a vector has index 1). The indices for events and stations are declared as arrays of integers, int idx_stat[N]; . Note that arrays are declared di erently than vectors; the size of the array is declared after the variable name. Constraints on variables are denoted by <lower=L, upper=U> a; , which means that the variable amust be larger than Land smaller than U. If a constraint is declared in the data block, Stan will check the values, and throw an error if the constraint is not met. Here, we declared constraints for the indices (they cannot be smaller than 1 or larger than the total number of events/stations), and the number of data/events/stations, since we assume that we will have at least one data point in the data. In theparameters {...} block, the parameters to be estimated are declared. As stated before, these are the coecients, the standard deviations, and the event and station terms. The parameters are declared similar as in the data block. Since the standard deviations are positive quantities, there must be a constraint <lower=0> in their declaration. The model model {...} block contains the probability model, which contains mainly the statements incrementing the log posterior. We have rst declared the prior distributions for the coecients and standard deviations, then for the event and station terms, and nally the likelihood. Note that the order of the statements is not important. Many functions in Stan are 9 vectorized, which means that they can be applied to vectors. The vectorization is a shorthand and more ecient than an explicit loop. Hence the statement c ~ normal(0,10); is equivalent to for(i in 1:K) { c[i] ~ normal(0, 10); } Similarly, the prior distributions for event/station terms, as well as the likelihood are vectorized. Next, we show an example how one could de ne the design matrix Xinside a Stan program. Here, we read in the predictor variables M,R, andVS30in thedata {...} block, and then de- ne the matrix Xin thetransformed data {...} block. The individual columns are assigned the linear predictors. Since most of the functions in Stan are vectorized, one can often assign the whole column at once. In the declaration of the linear prdictor accounting for magnitude dependent geometrical spreading, the term .*is used to indicate elementwise multiplication of two vectors. The rest of the model is unchanged. data { ... vector[N_eq] M; vector[N] R; vector[N_stat] VS; ... } transformed data { int K = 7; real h = 6; matrix[N,K] X; X[:,1] = rep_vector(1, N); X[:,2] = M[eq]; X[:,3] = square(8 - M[eq]); 10 X[:,4] = log(sqrt(square(R) + square(h))); X[:,5] = M[eq] .* log(sqrt(square(R) + square(h))); X[:,6] = R; X[:,7] = log(VS[stat] / 400); } parameters { ... } model { ... } 4.1 Running the Stan Model There are many di erent ways to execute the Stan model, depending on the working envi- ronment. Stan can be run from the command line, or from computing environments like R, Python, Julia, Matlab, Mathematica via dedicated packages (see https://mc-stan.org/ users/interfaces/ ). ForR(R Core Team, 2020), one can use the package Rstan (Stan De- velopment Team, 2020), or cmdstanR (Gabry and Cesnovar, 2020), which is a wrapper for the command line version. Similarly, for Python there exist the packages PyStan (Riddell et al., 2021) or cmdstanpy (Stan Development Team, 2021a). Rstan andPystan are more deeply integrated into Rand Python. We will show how to run the model from RusingcmdstanR . The reason for cmdstanR is that it can use the latest version of Stan, with the latest improvements. By contrast, Rstan lags behind in the versions of Stan used, but o ers more advanced features of analysis in R. In this section, we provide example code to run and assess a simple GMM using Stan. The full example is found at https://github.com/pstafford/StanGMMTutorial/blob/master/ gmm_stan_tutorial.html . Code to run the model using cmdStanR is shown below. data_list <- list(N = length(data_y), 11 N_eq = neq, N_stat = nstat, K = 7, X = data_x, Y = data_y, idx_eq = eqid, idx_stat = statid ) # sample file <- file.path( path/to/stan-file/ , gmm.stan ) mod <- cmdstan_model(file) fit <- mod$sample( data = data_list, seed = 123, chains = 4, iter_sampling = 1000, iter_warmup = 1000, refresh = 500 ) The Stan code is stored in a le called gmm.stan , and the command mod <- cmdstan_model(file) compiles the le. Then, with mod$sample() the MCMC sampler is called. In this case, we have speci ed a seed value for reproducibility. The number of chains is set to 4, with 1000 warm-up samples to be discarded, and 1000 samples post warm-up that can be used for in- ference. Every 500 iterations progress is printed. For an exhaustive list of arguments, see https://mc-stan.org/cmdstanr/reference/model-method-sample.html . The outcome of the model t is stored as a CmdStanMCMC object. The posterior package (B urkner et al., 2021) provides parctical tools to wok with the outcomes of Bayesian models, such as the draws format. The bayesplot package (Gabry et al., 2019) provides many plotting 12 phi_S2Sphi_SS tau 0200400600800100002004006008001000 020040060080010000.40.50.60.7 0.490.500.510.520.53 0.300.350.400.450.500.55Chain 1 2 3 4Figure 1: Traceplots of the posterior distributions of the standard deviation parameters. functions to assess the ts. For example, one can make a traceplot of the di erent chains for the standard deviation parameters using the following command posterior <- fit$draws() mcmc_trace(posterior, regex_pars = c("phi_SS", "tau", "phi_S2S"), facet_args = list(nrow = 2)) + xaxis_title(size = 30, family = "sans") + yaxis_title(size = 30, family = "sans") which produces a plot like in Figure 1. Histograms of the posterior distribution can be plotted with the function mcmc_hist() . An example is shown in Figure 2, which shows the posterior distributions of the coecients 13 c1;:::;c 7. This plot is created via mcmc_hist(posterior, regex_pars = c("c")) The plots from bayesplot are ggplot objects, and can be customized using functions from theggplot2 package (Wickham, 2016). In Figure 3, we show the histograms of the posterior distributions of the standard deviation parameters, together with the true values and the values estimated by lmer , using the following code df_lmer_sd <- as.data.frame(VarCorr(fit_lmer)) sd_names <- c( tau , phi_S2S , phi_SS ) size_text <- 20 size_title <- 30 i <- 1 # plot histogram for tau mcmc_dens(posterior, pars = sd_names[i]) + vline_at(c(sd_list[i], df_lmer_sd$sdcor[i]), size = c(1,0.75), linetype = c(1,2)) + xaxis_text(size = size_text, family = "sans") + yaxis_text(size = size_text, family = "sans") + xaxis_title(size = size_title, family = "sans") + yaxis_title(size = size_title, family = "sans") + grid_lines(color = "gray60") Athttps://github.com/pstafford/StanGMMTutorial , we provide some additional plots as- sessing the posterior distribution, such as pair plots and interval plots. For more possible options and examples, see https://mc-stan.org/bayesplot/index.html . We can also save the tted object with filename <- fit_gmm_stan.RDS fit$save_object(file = filename) Athttps://github.com/pstafford/StanGMMTutorial , we provide more on how one can re- trieve subsets of the parameters from the estimated posterior (in the form of a draws object), which can then be used for further analysis. 14 c[7]c[4] c[5] c[6]c[1] c[2] c[3] −1.0−0.50.0−3.50−3.25−3.00−2.75 0.240.280.320.36 −0.009−0.007−0.0055.07.510.012.515.017.5 −1.5−1.0−0.50.0−0.4−0.3−0.2−0.1Figure 2: Posterior histograms of the coecients. 01020304050 0.490.500.510.520.53 phi_SS 02468 0.4 0.5 0.6 0.7 tau 0.02.55.07.510.0 0.300.350.400.450.500.55 phi_S2S Figure 3: Posterior densities of standard deviation parameters, together with true values (solid vertical line) and values estimated by lmer (dashed vertical line) 15 4.2 Other Models Athttps://github.com/pstafford/StanGMMTutorial , we have compiled a large list of Stan models that showcase the rich variety of models that can be t using Stan. We also pro- vide the code to generate the data, run the model, and generate the Figures in this article. Models that can be t easily with Stan include truncated regression models (Kuehn et al., 2020b), models with correlated target variables (Kuehn and Scherbaum, 2015), measurement uncertainty models (Kuehn and Abrahamson, 2018; Sta ord, 2014), and models including spatial correlation Huang and Galasso (2019); Jayaram and Baker (2009); Kuehn and Abra- hamson (2020); Sta ord et al. (2019). Examples of these kinds of models can be found at https://github.com/pstafford/StanGMMTutorial , together with more information about implementation. In addition to these models, Kuehn (2021) provides code to estimate non- ergodic GMMs based on varying coecient models (Bussas et al., 2017; Gelfand et al., 2003; Landwehr et al., 2016). 5 Prior Distributions In our example, we have not put much emphasis on the prior distributions for the parameters. This is in part because we simulated enough data and the functional form is simple, so that the parameters are well constrained by the simulated data. In the development of GMMs, however, prior distributions can become important, as often data is sparse for some data ranges (e.g. at large magnitudes or very close to the fault), and in this case one needs to constrain the model behavior to avoid unphysical behavior. In GMMs estimated with maximum likelihood, some coecients are xed if they cannot be reliably estimated from data; however, this ignores possible uncertainty in the xed parameter. In terms of a Bayesian model, one can think of a xed parameter as a parameter with a very narrow prior distribution. It is dicult to give general advice about setting prior distributions that work in every case, but we believe that this is a par of the model that should no be neglected. Some general advice on prior distributions can be found at https://github.com/stan-dev/stan/wiki/ Prior-Choice-Recommendations . To set informative prior distributions, one can use for example models that were developed 16 for a di erent region (but a similar tectonic environment), making the assumption that ground- motion scaling for the same tectonic envionment is broadly similar. This approach has been taken in e.g. Kuehn and Scherbaum (2016) and Kowsari et al. (2020), where it was shown that the use of informative prior distributions stabilized the regression results for regions with sparse data. Another approach would be to use simulations to inform the prior distributions. Simulations have been used to constrain model parts that are not well informed by data, such as hanging wall e ects in the NGA West-2 models (Abrahamson et al., 2014; Boore et al., 2014; Campbell and Bozorgnia, 2014; Chiou and Youngs, 2014), based on Donahue and Abrahamson (2014). Similarly, the nonlinear site ampli cation in these models is based on simulations as well (Kamai et al., 2014; Walling et al., 2008). A similar approach was taken to constrain the large magnitude scaling in the NGA Subduc- tion model of Kuehn et al. (2020a), which is based on the models of Gregor et al. (2002) and Atkinson and Macias (2009). This is not well constrained by data, given that there are very few large magnitude events in the NGA Subduction database (Bozorgnia and Stewart, 2020), and thus it needs to be constrained based on external information. The near-fault term in the NGA subduction model of (Parker et al., 2020) is also based on simulations, though their model is not developed via Bayesian inference. Ordaz et al. (1994) derive prior distributions for a GMM applicable in Mexico City from seismological theory, and discuss in detail how they arrive at the values for the expectation and range of their prior distributions. In our experience, prior distributions are important for the standard deviations of par- tially nonergodic adjustment terms. Typically, a partially nonergodic model comprises few regions, so the standard deviations of the adjustment terms are not well constrained. In this case, it makes sense to use a prior distribution that discourages large values of the standard deviation, so that only a strong regional signal in the data leads to signi cant deviations of the regional parameters from the global average model. Kuehn et al. (2020a) compared Half- Cauchy and exponential prior distributions for ther regional standard deviations. They found that the Cauchy distributions, which has heavy tails, can lead to spurious regional devia- tions for some coecients, and concluded that the exponential distribution is a better choice. Inhttps://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations , sometimes 17 also a Half-normal or a Half-Student-T distribution with low degrees of freedom is recommended for standard deviations of hierarchical models. The notion discussed in the previous paragraph is informed by Simpson et al. (2017), who develop the penalized complexity (PC) prior approach, which sees a hierarchical model (such as a regional random e ects model) as a complex extension of a simpler base model. The more complex model is penalized (shrunk towards the base model) unless there is strong information in the data that supports it. This approach has been extended to priors for spatial models (Fuglstad et al., 2019), which becomes important for nonergodic models based on varying coecient models (Franco-Villoria et al., 2019; Landwehr et al., 2016). In the end, as stated before it is dicult to give general advice, but we believe that one should at least think about the prior distributions for each parameter. To get an idea about the impact of prior distributions, one can also simulate data based on the prior distributions (i.e. sample coecients from the prior distribution, then generate synthetic data based on the model). One can then assess whether the simulated data \makes sense", or show unrealistic values. This is a example of (qualitative) prior predictive checking (Gabry et al., 2019; Gelman et al., 2020). 6 Discussion We have provided a simple introduction to the estimation of GMMs with Stan, using MCMC. The main contribution, though, is the list of models compiled at https://github.com/pstafford/ StanGMMTutorial , which can serve as a starting point for anyone who is interested to esimate GMMs with Bayesian inference. Throughout the article, we have mentioned that we believe the Bayesian approach has advantages in GMM development: one can set informative prior distributions, and the posterior distribution allows for an easy quanti cation of epistemic un- certainty. There are disadvantages as well: a Stan model estimated with MCMC will generally take longer to run than an ML model. For example, the simple example model took about 5 minutes with Stan, but only about a second using lmer . For more complicated models, the runtime will typically be even longer. It depends on the actual application and model whether this is a concern. This article is only a brief introduction to Stan, and touches upon Bayesian statistics very 18 lightly. For more thorough exposition, some introductory textbooks, such as Kruschke (2015) or McElreath (2020) are recommended. For a better understanding of Stan, one should also look at the Stan documentation ( https://mc-stan.org/users/documentation/ ). There are di erent ways to code the same model in Stan, and it is often not immediately obvious which pa- rameterization is most ecient. An example is the di erence between a centered/non-centered parameterization in hierarchical models (Betancourt, 2020; Betancourt and Girolami, 2015). For model comparison in MCMC models, one can use the widely applicable information criterion (WAIC) (Watanabe, 2013), cross-validation (CV), or approximate leave-one-out CV (LOO-CV) (Vehtari et al., 2017). Many of these are easily available for Stan models, and can be computed via the looR-package (Vehtari et al., 2020a). We have focused on running the Stan models in R, but Stan can be run from many interfaces. For Python users, there exist the packages PyStan orcmdstanpy . For assessment of the posterior dstribution, the Python package ArviZ (Kumar et al., 2019) is very useful, and provides many of the same tools and capabilities as the R-packages used in this work.