## 5 Discussion  
This paper proposes a novel hierarchical GLM framework to estimate insurance risk premium with the incorporation of the heterogeneity and NDB covariate, outlined in research questions RQ1 and RQ2. Our proposed model presents overall good empirical performance results in retrieving the original error-free parameter values, which suggests that it is worth considering the Gustafson correction boosted with the hierarchical Bayesian framework to avoid the impact of the NDB errors. Ultimately, we seek to develop practical guidelines for selecting the optimal scaling factor with the Gustafson correction in the absence of gold standard data.  
Regarding **RQ1**, we initially analyzed the partial pooling effect of a hierarchical GLM as a solution for addressing heterogeneity. Given the predetermined risk clusters \(j=1,\ldots,6\), we observed the highest LPPD value from the gold standard model, indicating an improvement in predictive power compared to other conventional GLM models. This suggests that the effect of heterogeneity is effectively mitigated. For **RQ2**, we expanded the model by incorporating Gustafson correction into the hierarchical GLM framework to mitigate the model risk linked to the NDB covariate. Throughout our experiments, we demonstrated the impact of the model risk when using a naive model built on the NDB covariate. This model misspecification led to significant distortion in the spread and created extreme modality in the predictive distribution. In this context, we showed how applying the Gustafson correction within the same framework effectively'restores' the predictive distribution's original spread and characteristics. The NDB error mitigation was further confirmed by additional performance metrics such as LPPD, SSPE, SAPE, \(D_{KL}\), etc.  
A key aspect of our hybrid modeling framework to handle NDB covariates is the choice of scaling factor \(\zeta\), which our findings demonstrate to be crucial in error correction.  
Figure 8: Curve alignment under the condition of ‘error rate = 0.40’ with \(\zeta=0.5\): cluster-wise histograms (for \(j=1,\ldots 6\)) of the observed claim amount \(Y_{h}\) on a log scale and the out-of-sample predictive densities obtained from Model (A), -(B), and -(C).  
Specifically, \(\zeta\) moderates the contribution of \(\widehat{\lambda}^{2}:V(\mathbf{x}^{*}|\mathbf{z})\) to explaining \(\mathbf{\mathfrak{T}}^{2}:V(\mathbf{x}^{*}|\mathbf{x})\), guiding adjustments of the erroneous parameters. Notably, we found that both the scaling factor \(\zeta\) and the error rates \(R_{e_{x}}\) in the NDB covariate are crucial, requiring careful calibration for optimal results. To determine the ideal \(\zeta\), we devised a rule of thumb relying on the maximum LPPD value, providing a path forward without access to gold-standard data. For the error rates \(R_{e_{x}}\) in the NDB covariate, our experimental results indicate that the hierarchical GLM equipped with the Gustafson correction performs consistently well across different error rate scenarios, particularly within the range of scaling factor such as \(0.5\leq\zeta\leq 0.7\). However, these numbers suggest a relatively strong degree of relation between \(\widehat{\lambda}^{2}:V(\mathbf{x}^{*}|\mathbf{z})\) and \(\mathbf{\mathfrak{T}}^{2}:V(\mathbf{x}^{*}|\mathbf{x})\), as reflected by the prior knowledge: \(\mathbf{\mathfrak{T}}_{j}^{2}=(1-\zeta)\widehat{\lambda}_{j}^{2}\) in Equation (34). Furthermore, the optimal range of \(\zeta\) appears to be dataset dependent, as it is likely unique to each one. Varying relationships among covariates can lead to distinct ranges for \(\zeta\). Indeed, we have yet to explore other datasets where this degree of relation is characterized by a different range of \(\zeta\). It can be important to investigate how the performance of the correction varies across different error rate scenarios when the optimal \(\zeta\) is ranged differently. This leaves room for further research to fully understand the implications of \(\zeta\) on model performance.