## 3 Modeling Method for \(\mathbf{S_{h}}|\mathbf{X}^{F},\mathbf{X}^{S}\)
### Introduction of Gustafson Correction to Address NDB Covariate: RQ2
Algorithm S2 in Part 2 of the Supplementary Materials provides a detailed description of the Gibbs sampler tailored for modeling claim amounts \(Y_{h}\) using the log-normal density. The approach for claim counts \(N_{h}\) based on the negative binomial distribution would be analogous and is therefore omitted for brevity. At the initialization of the Gibbs sampler, it is crucial to establish the initial values for the hyperparameters \(\mathbf{\phi}:\{\mathbf{\beta}_{0}^{S},\Sigma_{\beta_{0}}^{S},u_{0}^{S},v_{0}^{S},m_{0}, \delta,q_{0},\Lambda,\rho_{u1},\rho_{u2},\rho_{v1},\rho_{v2},\epsilon_{0}^{S},d_{0 }^{S},S_{0}^{S},h_{0}^{S}\}\), which support the prior choices in Equations (14c) and (14e) for the log-normal GLM. Using these values, the initial parameters for the outcome \(\mathbf{\theta}^{(old)}:\{\mathbf{\beta}_{j}^{S(old)},\sigma_{j}^{2(old)}\}\) and covariate \(\mathbf{w}:\{\mathbf{\tau}_{j}^{S},\lambda_{j}^{2S}\}\) are determined. These parameters ultimately yield the values for the communal hyperparameters \(\mathbf{\beta}_{0}^{S+},\Sigma_{\beta_{0}}^{S+},u_{0}^{S+},v_{0}^{S+}\), which guide both complete pooling and no-pooling throughout the Gibbs sampling process. The two-stage Gibbs sampler for the log-normal model is executed as follows:  
**[Stage.1] Sampling with Complete Pooling**  
The Gibbs sampler first estimates the outcome parameters \(\mathbf{\theta}\) (without clustering) using a complete pooling approach. This aims to refine the communal hyperparameters \(\mathbf{\beta}_{0}^{S+},\Sigma_{\beta_{0}}^{S+},u_{0}^{S+},v_{0}^{S+}\) as outlined in Equations (23) and (24). Given the absence of conjugate priors for the log-normal outcome, the Metropolis-Hastings (MH) algorithm is used,  
\[\beta_{0}^{F}\mid\underline{m_{0}},\underline{\measuredangle},N,\mathbf{X}^{F}, \ \Sigma_{\beta_{0}}^{F}\sim\ \mathbf{MVN}\Bigg{(}\frac{\hat{\xi}}{\hat{\xi}+1} \underline{m_{0}}+\frac{1}{\hat{\xi}+1}\beta^{F},\ \frac{\Sigma_{\beta_{0}}^{F}}{\hat{\xi}+1}\Bigg{)}\]1 The acyclic graphical representation of the flows of the parameter updates in the hierarchical GLM. This is a snapshot for a single iteration (M = 1). with the prior serving as its proposal distribution, from where the candidate values are drawn. Proposal samples \(\mathbf{\theta}^{(new)}\) are generated based on these communal hyperparameters, yielding \(\mathbf{\theta}^{(*)}\), which support the global mean \(E[Y_{hi}|\mathbf{X}]\), while re-updating the communal hyperparameters \(\mathbf{\beta}_{0}^{S+},\Sigma_{\beta_{0}}^{S+},u_{0}^{S+},v_{0}^{S+}\). See Algorithm S2 in Part 2 of the Supplementary Materials. **[Stage.2] Sampling without Pooling**  
In this stage, with the cluster membership \(j\) already established, the goal is to generate accurate parameter estimates \(\mathbf{\theta}^{(*)}\) for each cluster. These estimates are guided by the communal hyperparameters \(\mathbf{\beta}_{0}^{S+},\Sigma_{\beta_{0}}^{S+},u_{0}^{S+},v_{0}^{S+}\) obtained from **[Stage.1]**. By using these communal hyperparameters, the Gibbs sampler optimizes \(\mathbf{\theta}_{j}^{(*)}\) for each risk cluster, minimizing the within-cluster variability. The log-likelihood value can be computed after each iteration to monitor convergence. The implementation detail is provided in Algorithm S2 in Part 2 of the Supplementary Materials. Clustering \(S_{h}|\mathbf{X}^{F},\mathbf{X}^{S}\) with NDB Case Covariate: RQ2 Tackling NDB Errors  
This section introduces a novel method for handling the NDB covariate (RQ2) based on the hierarchical GLM framework and partial pooling for risk premium modeling. Consistent with the parametric Bayesian principle [28], we assume that the risk clusters \(j=1,\ldots J\) have already been established.