Dirichlet Process Log Skew-normal Mixture with Missing at Random Covariate in Insurance Claim Analysis

1. Introduction

In short-term insurance contracts, predicting accurate aggregate claims is essential for major actuarial decisions such as pricing or reserving. However, it is often not easy to model the aggregate loss due to its complex distributional features such as high skewness, zero inflation, hump shape, multi-modality, etc. With the advance of the modern Bayesian paradigm and computing power, the development of full distribution of aggregate claims has been studied and applied in actuarial practice. In particular, because of its considerable flexibility, a Bayesian nonparametric (BNP) approach has been gradually recognized to solve distributional conundrums in an insurance context. For instance, Hong and Martin (2018) [1] recently developed the Dirichlet process model as a BNP approach that maximizes the fitting flexibility of the full distribution for insurance loss, which obviates the chance of model misspecification bias. In this article, as an extension of their work, we aim to go beyond the search for the maximized fitting flexibility, focusing on the issues that arise from the presence of covariates and the aggregate outcome (total losses). The implication is that the predictive distribution for the expected aggregate claims developed under Hong and Martin’s Dirichlet process framework can be biased with the incorporation of covariate effects and log-normal convolution. For example, as covariates add new information that differentiates the data points of the outcome variable, a new structure can be introduced into the data space, and this increases the within-cluster heterogeneity [2]. Besides, the incorporation of missing covariates may exacerbate the existing heterogeneity. Additionally, assuming that the outcome variable describes the aggregate losses, rather than individual claim amounts, it is difficult to compute the log-normal convolution as it does not have a closed-form solution. In this regard, our study extends their work by addressing the following research questions:

• RQ1. If an additional unobservable heterogeneity is introduced by the inclusion of covariates, what is the best method to capture the within-cluster heterogeneity in modeling the total losses, comparing several conventional approaches?
• RQ2. If an additional estimation bias results from the use of the incomplete covariates under Missing At Random (MAR), what is the best way to increase the imputation efficiency, comparing several conventional approaches?
• RQ3. If an individual loss is distributed with log-normal densities, what is the best way to approximate the sum of log-normal outcome variables, comparing several conventional approaches?

2. Discussion on Research Questions and Related Work

Let $Y_i,i=1,2,\dots ,N$ be the independent claim amount (reported by each policyholder for a single policy) random variable, defined on a common probability space ($\mathrm{\Omega },\mathcal{F},P$) from a certain loss distribution such as log-normal. Let $\mathit{X}$ be a vector of covariates, and $N\left(t\right)$ be the total claim count denoting the number of individual claims for a single policy up to time t (policy period). The aggregate claim $S_h\left(t\right)$ for a single policy, h, given time t can be expressed as a convolution: $S_h\left(t\right)={\sum }_{i=1}^{N\left(t\right)}{Y}_i=Y_1+Y_2+\dots +Y_{N\left(t\right)}$. At the end of the policy period t, let $\tilde{S}\left(t\right)$ be the total aggregate claim amounts from the total policies received by an insurer, then: $\tilde{S}\left(t\right)={\sum }_{h=1}^H{S}_h\left(t\right)=S_1\left(t\right)+S_2\left(t\right)+\dots +S_H\left(t\right)$ in which H is the total number of policies on the contracts. Note that both convolutions described so far are built upon the assumption that the summands - $Y_i,i=1,2,\cdots ,N\left(t\right)$ and $S_h,h=1,2,\cdots ,H$ - are mutually independent and identically distributed with log-normal densities (to maintain homogeneity of each loss).

However, the involvement of covariates and the lack of closed-form solutions for the log-normal sum bring about several challenges that violate the assumptions for an accurate estimation of the total aggregate loss $\tilde{S}\left(t\right)$. To begin with, the use of covariates gives rise to an additional within-cluster heterogeneity. Kass et al.(2008) [3] describes a standard aggregate loss modeling principle denoting that the expected aggregate claims $E\left[S_h\right]$ is obtained by the product of the mean claim counts and severities: $E\left[S_h\right]=E\left[N\right]E\left[Y\right]$. With the inclusion of covariates $\mathit{X}$, a new unknown structure or heterogeneity is introduced into the data space of $Y_i$, which means that $Y_1|{\mathit{X}}_1,Y_2|{\mathit{X}}_2,\cdots ,Y_N|{\mathit{X}}_N$ within a single policy can still be independent, but cannot be identically distributed. Therefore, $E\left[S_h\right|\mathit{X}]\ne E\left[N\right|\mathit{X}]E\left[Y\right|\mathit{X}]$, and the total aggregate loss $\tilde{S}\left(t\right)$ becomes difficult to compute with the conventional collective risk modeling approach. In addition, assuming that the severity $Y_i$ follows a log-normal distribution, the computation of $\tilde{S}\left(t\right)$ becomes quite difficult as its convolution $S_h$ is not known to have a closed-form [4]. Another challenge is the missing covariates in $S_h|\mathit{X}$. As shown by Ungolo et al.(2020) [5], the missing covariates under the missingness at random (MAR) assumption lead to the biased parameter estimations because the uncertainty in the estimation results of the parameters describing the outcome Y is heavily affected by the quality of covariates $\mathit{X}$. Again, in this case, $\tilde{S}\left(t\right)$ cannot be computed properly.

Compounding all this, we propose the Dirichlet process log skew-normal mixture to model the $S_h|\mathit{X}$. We consider the Dirichlet process framework to cope with the within-cluster heterogeneity as suggested by Hong and Martin (2018); Braun et al.(2016) [1,6] while employing the log skew-normal approximation studied by Li (2008) [7] to compute each $S_h|\mathit{X}$, the sum of log-normal random variables ${\sum }_{i=1}^{N\left(t\right)}{Y}_i|\mathit{X}$. When it comes to the problem of missing covariates, we exploit the generative capability of the Dirichlet process to capture the latent structure of data, which allows for a rigorous statistical treatment of MAR covariates.

2.1. Can Dirichlet Process Capture the Heterogeneity and Bias?: RQ1, RQ2

In Figure 1, $Y_i$ refers to individual claim amount and each $S_h$ represents a total claim amount defined by a unique policy (cluster h) as a homogeneous distribution. Although an insurer can collect the aggregate loss data $S_h$ for each policy cluster given policy period t, individual policyholders (in different risk classes) can raise more than one claim (i.e. random $N_h\left(t\right)$) at any time over a fixed time horizon t, and their corresponding claim amounts (i.e. random $Y_h^{*}$) will not be known in advance. Hence, the unsettled liability information of $Y_h^{*}$ from certain policyholders always renders $S_h$ incomplete, which is often translated into the challenge of their inherent stochasticity. In addition, the new claims $Y_h^{*}$ raised from unknown risk classes can trigger inherent heterogeneity across unique clusters as well. To make matters worse, if introducing covariates $\mathit{X}$ to better understand the different risk classes, one might introduce an additional source of heterogeneity into the scene, which prevents each cluster from being identically distributed.

With respect to this, Hong and Martin (2018) [1] propose the concept of the loss distribution mixture for each cluster based on the Dirichlet process framework. The main idea behind the Dirichlet process mixture (DPM) is to produce a single master distribution to model stochasticity in $S_h$ with the help of an infinite dimensional parametric structure and the probabilistic simulations of clustering scenarios. Braun et al.(2006) [6] articulates how the DPM automatically captures unobservable heterogeneity such as intracorrelation between claim amounts $Y_i$ in the different risk classes without specifying the number of the classes upfront. In short, no matter how complex the distribution of the data is, the DPM is capable of accommodating any distributional properties - multi-modes, skewness, heavy tails, etc. - resulting from unobservable heterogeneity; and therefore, dramatically minimizes model misspecification biases.

With the inclusion of the covariates, the DPM offers a useful bedrock for a MAR treatment. As a generative modeling approach, the DPM models both outcomes $S_h$ and covariates $\mathit{X}$ jointly to produce cluster memberships. This is used as key knowledge to identify the latent structure of the data. For example, in the domain of medicine research, Roy et al. (2018) [8] develop a novel imputation strategy for the MAR covariate, using the latent structure unraveled by the DPM and the other covariate knowledge available. A further survey of imputation methods based on the Nonparametric Bayesian framework can be found in Si and Reiter (2013) [9] and references therein.

2.2. Can Log Skew-Normal Mixture Approximate the Log-Normal Convolution?: RQ3

The log-normal distribution has been considered a suitable claim amount $Y_i$ distribution due to its non-negative support, right-skewed curve, and moderately heavy tail to accommodate some outliers. However, if generalizing the individual claim amount $Y_i$ by introducing a log-normal distribution, the convolution computation for $S_h$ fails because the exact closed form for the log-normal sum is unknown.

Furman et al. (2020) [10] present several existing methods for the log-normal sum approximation that have been studied in the literature. This includes the moment matching approximation approaches such as Minimax approximation, Least squares approximation, Log shifted gamma approximation, and Log skew-normal approximation. The distance minimization approaches - Minimax approximation or Least squares approximation - described by Beaulieu and Xie (2003); Zhao and Ding (2007) [4,11] are conceptually simple, but they require to fit the entire cumulative densities to the sum of claim amounts, which can be computationally expensive and easy to fail when the number of the summands $Y_i$ increases. The Log shifted gamma approximation suggested by Lam and Le-Ngoc (2007) [12] has less strict distributional assumptions, but it is not very accurate at the lower region of the distribution. In our study, special attention is paid to the possibility of the Log skew-normal approximation method for the sake of simplicity. A skew-normal distribution as an extension of a normal distribution has a third parameter to naturally explain skewness apart from the other parameters (for a location and spread). Li (2008) [7] points out that one can exploit the third parameter of the skew-normal distribution to capture different skewness levels of each summand. Taking the log of skew-normal densities, we can approximate $S_h$, the sum of the log-normal $Y_i$. Using the log skew-normal as the underlying distribution for $S_h$ in the DPM framework, one can eliminate the need to compute the cumulative density curve, and its closed-form density and the optimal distribution parameters for $S_h$ can be easily obtained by the moment matching technique. For further details, see Li (2008) [7] and the references contained within.

2.3. Our Contribution and Paper Outline

The contribution of this study is as follows: first, using the Bayesian nonparametric framework, we propose solutions to the two major challenges of the aggregate claim $S_h$ computation - 1) heterogeneity in the log-normal random variable $Y_i$, 2) lack of closed-form of the sum of log-normal random variables $Y_i$ - in a more unified fashion. Second, we introduce covariates $\mathit{X}$ into the aggregate claim modeling framework, taking into account the adverse impact triggered by the covariates $\mathit{X}$. This includes the added heterogeneity across $Y_i$ and the missing information fed by MAR covariates $\mathit{X}$. To our knowledge, there have been no previous attempts to either estimate the log skew-normal mixture within the DPM framework or use the DPM to handle the MAR covariate in the insurance loss modeling.

The rest of the paper is structured as follows. In Section 3, we describe the proposed modeling framework for $S_h$, assuming log-normal distributed $Y_i$ and the inclusion of both continuous and discrete covariates $\mathit{X}$. This section also presents our novel imputation approach for the MAR covariate within the DPM framework. Section 4 clarifies the final forms of the posterior and predictive densities accordingly. Section 5 presents our empirical results, and validates our approach by fitting to two different datasets with different sample sizes drawn from the R package CASdatasets and the Wisconsin Local Government Property Insurance Fund (LGPIF). This is followed by a discussion in Section 6.

3. Model: DP Log Skew-Normal Mixture for $S_h|X$

3.1. Background

Consider that there are multiple unknown risk classes (clusters) across the claim $Y_i$ information within each policy, and then the individual aggregate claims $S_h$ for the policy h would have diverse characteristics that cannot be explained by fitting a single log skew-normal distribution. In order to approximate the distribution that captures such diverse characteristics in $S_h$, we seek to investigate diverse clustering scenarios. To this end, as suggested by Hong and Martin (2018) [1], we exploit the infinite mixture of log skew-normal clusters and their complex dependencies by employing a Dirichlet process. The Dirichlet process produces a distribution over clustering scenarios (with clustering parameters).

$$\begin{array}{cc}\hfill \{{\mathit{\theta }}_j,{\mathit{\omega }}_j\}& \sim G\hfill \\ \hfill G& \sim DP\left(\alpha ,G_0\right)\hfill \end{array}$$

where G denotes the clustering scenarios, and the important components of G are

• ${\mathit{\theta }}_j$: the parameters of the outcome variable defined with cluster j.
• ${\mathit{\omega }}_j$: the parameter of the cluster weights defined with cluster j.

G, as a single realization of the joint cluster probability vector $\{G\left(A_1\right),G\left(A_2\right)\dots \}$ sampled from the DPM model, takes independent partitions $A_1,A_2,\dots$ of the sample space ${\bigcup }_{k=1}^{\infty }{A}_k=A$ of the support of $G_0$. By sufficient simulations of G, the Dirichlet process investigates all possible clustering scenarios rather than relying on a single best guess. The overall production of G is controlled with two parameters - a precision $\alpha$ and a base measure $G_0$. The precision $\alpha$ controls a variance of sampling G in the sense that larger $\alpha$ generates new clusters more often to account for the unknown risk classes. The base measure $G_0$, as the mean of DP($\alpha ,G_0$), is a DP prior over the joint space of all parameters for the outcome model, covariate model, and the precision $\alpha$, as shown in Ghosal (2010) [13].

Note that the original research on DPM by Hong and Martin (2018) [1] mainly focuses on the random cluster weights ${\mathit{\omega }}_j$ independent of the covariates $\mathit{X}$. On the other hand, in our model, the covariate effects are incorporated into the development of cluster weights ${\mathit{\omega }}_j$. All calculations for the development of the DPM modeling components in this paper are based on the principles introduced by Ferguson (1973), Antoniak (1974), and Sethuraman (1994) [14,15,16].

3.2. Model Formulation with Discrete and Continuous Clusters

Let the outcome be $S=\{S_1,S_2,\dots ,S_H\}$ denoting the H different aggregate claims (incurred by the H different policies). We assume that the covariate ${\mathit{x}}_1$ is binary, and the ${\mathit{x}}_2$ is Gaussian, and then our baseline DPM model can be expressed as:

$$\begin{array}{cc}\hfill S_h|{\mathit{x}}_1,& {\mathit{x}}_2,{\mathit{\beta }}_j,{\sigma }_j^2,{\xi }_j,{\tilde{\mathit{\beta }}}_j\hfill \\ \hfill & \sim \delta \left({\mathit{X}}^T{\tilde{\mathit{\beta }}}_j\right)\mathbb{1}(S_h=0)+\left[1-\delta \left({\mathit{X}}^T{\tilde{\mathit{\beta }}}_j\right)\right]LogSN\left({\mathit{X}}^T{\mathit{\beta }}_j,{\sigma }_j^2,{\xi }_j\right)\hfill \\ \hfill {\mathit{x}}_1|{\pi }_j& \sim Bern\left({\pi }_j\right)\hfill \\ \hfill {\mathit{x}}_2|{\mu }_j,{\tau }_j^2& \sim N\left({\mu }_j,{\tau }_j^2\right)\hfill \\ \hfill \{{\mathit{\theta }}_j,{\mathit{w}}_j\}& \sim G\hfill \\ \hfill G& \sim DP\left(\alpha ,G_0\right)\hfill \end{array}$$

(1)

where j is the risk class index; ${\mathit{\theta }}_j=\{{\mathit{\beta }}_j,{\sigma }_j^2,{\xi }_j,{\tilde{\mathit{\beta }}}_j$} describe the outcome model while ${\mathit{w}}_j=\{{\pi }_j,{\mu }_j,{\tau }_j^2\}$ explains the covariate model. $S_h$ is modeled as a mixture of a point mass at 0 and positive values distributed with log skew-normal density to address the complications of zero inflation in the loss data. $\delta \left({\mathit{X}}^T{\tilde{\mathit{\beta }}}_j\right)$ models the probability of the outcome being zero using a multivariate logistic regression. Variable Definitions section has a brief description of all parameters used in this study.

Considering a Dirichlet process log skew-normal mixture to house the multiple unknown risk classes in $S_h$, it is necessary to differentiate the forms of mixture components depending on the types of clusters it uses - the discrete and continuous. While keeping the inference of the cluster parameters to be data dominated, the DPM first develops discrete clusters based on the given claim information and then extrapolates certain unobservable clusters of claims by examining the heterogeneity (or hidden risk classes) of each cluster. In this process, the DPM develops new continuous clusters additionally and assesses them with some probabilistic decision-making algorithms, rendering the parameter estimations computationally efficient and asymptotically consistent [17].

The discrete mixture components (clusters) in the DPM framework have the standard form that is useful in accounting for the observed classes such as policy information for aggregate loss $S_h$ [18]. In calculating the discrete cluster probabilities, we assume that the non-zero outcome and covariates are distributed with the densities denoted by

$$\begin{array}{cc}\hfill f_{LSN}\left(S_h|{\mathit{X}}_h^T{\mathit{\beta }}_j,{\sigma }_j^2,{\xi }_j\right)& =\frac{2}{S_h{\sigma }_j}\varphi \left(\frac{log{S}_h-{\mathit{X}}_h^T{\mathit{\beta }}_j}{{\sigma }_j}\right)·\mathrm{\Phi }\left({\xi }_j·\frac{log{S}_h-{\mathit{X}}_h^T{\mathit{\beta }}_j}{{\sigma }_j}\right)\hfill \end{array}$$

(2a)

$$\begin{array}{cc}\hfill f_{Bern}\left({\mathit{x}}_1|{\pi }_j\right)& ={\pi }_j^{{\mathit{x}}_1}{\left(1-{\pi }_j\right)}^{1-{\mathit{x}}_1}\hfill \end{array}$$

(2b)

$$\begin{array}{cc}\hfill f_N\left({\mathit{x}}_2|{\mu }_j,{\tau }_j^2\right)& ={\displaystyle \frac{1}{\sqrt{2\pi {\tau }_j^2}}}exp\left\{-{\displaystyle \frac{1}{2{\tau }_j^2}}{\left({\mathit{x}}_2-{\mu }_j\right)}^2\right\}\hfill \end{array}$$

(2c)

where $\varphi (·)$ and $\mathrm{\Phi }(·)$ are standard normal probability and cumulative density functions for the log skew-normal density. To model the outcome data $S_h|{\mathit{X}}_h$ for the policy h, the DPM takes the general form of the mixture

$$f\left(S_h\right|{\mathit{X}}_h,\mathit{\theta })=\sum _{j=1}^J{\mathit{\omega }}_j\left(\delta \left({\mathit{X}}_h^T{\tilde{\mathit{\beta }}}_j\right)\mathbb{1}(S_h=0)+\left[1-\delta \left({\mathit{X}}_h^T{\tilde{\mathit{\beta }}}_j\right)\right]f_{LSN}\left(S_h\right|{\mathit{X}}_h,{\mathit{\theta }}_j)\right)$$

(3)

where J is the total number of mixture components (risk classes), ${\mathit{\theta }}_j=\{{\mathit{\beta }}_j,{\sigma }_j^2,{\xi }_j,{\tilde{\mathit{\beta }}}_j\}$ and ${\mathit{w}}_j=\{{\pi }_j,{\mu }_j,{\tau }_j^2\}$ are the outcome and covariate parameters to explain the risk clusters, and ${\mathit{\omega }}_j$, functions of covariates: ${\mathit{\omega }}_j\left({\mathit{X}}_h\right|{\mathit{w}}_j)$, are the cluster components weights (mixing coefficient) satisfying ${\sum }_{j=1}^J{\mathit{\omega }}_j=1$.

However, when the DPM is extended as $j\to \infty$, the new continuous clusters are introduced by the $G_0$ (with its infinite-dimensional parametric structure) in order to address the additional unknown risk classes. This assesses the within-class heterogeneity in $S_h$ by confronting the current discrete clustering result and investigating the homogeneity more closely. As the new clusters are considered countably infinite, their corresponding forms of the outcome and covariate models to obtain the continuous cluster are given by

$$\begin{array}{cc}\hfill f_0\left(S_h\right|{\mathit{X}}_h)& =\int f\left(S_h\right|{\mathit{X}}_h,\mathit{\theta })d{G}_0\left(\mathit{\theta }\right)\hfill \end{array}$$

(4a)

$$\begin{array}{cc}\hfill f_0\left({\mathit{x}}_1\right)& =\int f_{Bern}\left({\mathit{x}}_1\right|\mathit{w})d{G}_0\left(\mathit{w}\right)\hfill \end{array}$$

(4b)

$$\begin{array}{cc}\hfill f_0\left({\mathit{x}}_2\right)& =\int f_N\left({\mathit{x}}_2\right|\mathit{w})d{G}_0\left(\mathit{w}\right)\hfill \end{array}$$

(4c)

They are also known as a “parameter-free outcome model" and a “parameter-free covariate model" respectively to develop the new continuous cluster mixture. Given a collection of outcome-covariate data pairs $D={\{S_h,{\mathit{X}}_h\}}_{h=1}^H$, the DPM puts together the current discrete clusters and new continuous clusters to update the mixture form in Equation (3), with help of Monte Carlo Markov Chain (using sufficiently simulated samples of the major parameters ${\mathit{\theta }}_j,{\mathit{w}}_j$). Consequently, the sample G described in Equation (1) becomes $G=f\left(S_h\right|{\mathit{X}}_h,D)={\sum }_{j=1}^{\infty }{\mathit{\omega }}_j·{\delta }_{z_j}$ where ${\delta }_{z_j}$ denotes both discrete and continuous cluster densities as point mass distributions at the random locations sampled from $G_0$. Aligned with such flexible cluster development, the form of the predictive distribution can be molded based on the knowledge extracted from G, as follow:

$$f\left(S_h\right|{\mathit{X}}_h,\mathit{\theta },\mathit{w},\alpha )=\frac{{\mathit{\omega }}_{J+1}^{*}}{{\mathit{\omega }}_{J+1}^{*}+{\sum }_{j=1}^J{\mathit{\omega }}_j^{*}}·f_0\left(S_h\right|{\mathit{X}}_h)+\frac{{\sum }_{j=1}^J{\mathit{\omega }}_j^{*}·f\left(S_h\right|{\mathit{X}}_h,{\mathit{\theta }}_j)}{{\mathit{\omega }}_{J+1}^{*}+{\sum }_{j=1}^J{\mathit{\omega }}_j^{*}}$$

(5)

and the finalized cluster weights in Equation (5) are secured through computing these two sub-models below for discrete and continuous cluster weights respectively which reflect the properties of the clusters and relevant covariates.

$$\begin{array}{cc}\hfill {\mathit{\omega }}_{J+1}^{*}& =\frac{\alpha }{\alpha +H}·f_0({\mathit{x}}_1,{\mathit{x}}_2)\hfill \end{array}$$

(6a)

$$\begin{array}{cc}\hfill {\mathit{\omega }}_j^{*}& =\frac{n_j}{\alpha +H}·f({\mathit{x}}_1,{\mathit{x}}_2|{\mathit{w}}_j=({\pi }_j,{\mu }_j,{\tau }_j^2))\hfill \end{array}$$

(6b)

where $\alpha$ is the precision parameter to control the acceptance chances of the new clusters, $n_j$ is the number of observations in cluster j, $f_0\left(\mathit{X}\right)$ is the parameter-free covariate model in Equation (4b, 4c) to support the new continuous clusters, and $f\left(\mathit{X}\right|{\mathit{w}}_j)$ is the covariate model to support the current discrete clusters. Note that instead of the popular stick-breaking scheme used by Hong and Martin (2018) [1], the cluster weights are obtained based on the covariate models of ${\mathit{x}}_1,{\mathit{x}}_2$ that explain the outcome $S_h$.

The simulated outcome model $f\left(S_h\right|{\mathit{X}}_h,D)={\sum }_{j=1}^{\infty }{\mathit{\omega }}_j·{\delta }_{z_j}$ and its predictive model in Equation (5) show that although the DPM framework allows infinite-dimensional modeling, the dimension of the sampling output G is adaptive as it is a mixture with at most finite components determined by data itself (e.g. its dimension cannot be greater than the total sample size H). This gives the model flexibility, and throughout such modeling flexibility, the G can become the comprehensive mixture distribution for $S_h$, accommodating all distributional properties of the given claims as well as the additional unknown claims.

3.3. Modelling $S_h|X_h$ with Complete Case Covariate

The joint posterior update for the outcome and covariate parameters - ${\mathit{\theta }}_j,{\mathit{w}}_j$ - in Equation (5,6) can be made through a Gibbs Sampler. Using the conditional distribution of the unobservable variables given the observed data, the Gibbs sampler can obtain draws from the analytically intractable posterior distribution of the parameters [20]. Let the cluster-index $j=1,2,\cdots ,J$ for the observation h be $s_h$. The parameter inference steps to ensure convergence are described below.

Step.1

Initialize the cluster membership and the main parameters:

(a)

First the cluster membership $j=1,\cdots ,J$ is initialized by some clustering methods such as hierarchical clustering or k-means, etc. This step provides an initial clustering of the data ($S_h,{\mathit{X}}_h)$ as well as the initial number of clusters.

(b)

Next, after all observations have been assigned to a particular cluster $j=1,2,\cdots ,J$, we can then update the parameters $\alpha$ and (${\mathit{\theta }}_j,{\mathit{w}}_j$) for each cluster. This is done using the posterior densities denoted by $p\left(\alpha \right|J),p\left(\mathit{\theta }\right|S_h,{\mathit{X}}_h)$, and $p\left(\mathit{w}\right|{\mathit{X}}_h)$ in which ($S_h,{\mathit{X}}_h$) represent all observations in cluster j.

Step.2

Loop through the Gibbs sampler and new continuous cluster selection:

Once the cluster memberships and parameters are initialized, we then loop through the Gibbs sampler many times (e.g. $M=100,000$ iterations) where the algorithm alternates between updating the cluster membership for each observation and updating the parameters given the cluster partitioning. Each iteration might give a slightly different selection of the new clusters based on the Polya Urn scheme [20], but the log-likelihood calculated at the end of each iteration can help keep track of the convergence of the selections. A detailed description of each iteration is given in Algorithm (A2) in Appendix B. The term $p\left(s_h\right|s_{-h})$ on lines 6 and 9 in Algorithm (A2) is the Chinese Restaurant process [19] posterior value given by

$$p\left(s_h\right|s_{-h})=\left\{\begin{array}{cc}c·{\displaystyle \frac{n_j^{-h}}{\alpha +H-1}},\hfill & \mathrm{for}\mathrm{record}h\mathrm{entering}\mathrm{into}\mathrm{existing}\mathrm{cluster}\text{:}{s}_h=j.\hfill \\ c·{\displaystyle \frac{\alpha }{\alpha +H-1}},\hfill & \mathrm{for}\mathrm{record}i\mathrm{entering}\mathrm{into}\mathrm{a}\mathrm{new}\mathrm{cluster}\text{:}{s}_h=J+1.\hfill \end{array}\right.$$

(7)

where c is a scaling constant to ensure that the probabilities sum to 1, and $s_{-h}$ is the collection of cluster indices $(s_1,s_2,\cdots ,s_{h-1},s_{h+1},\cdots ,s_H)$ assigned to every observation without the cluster index $s_h$ of the observation h. As shown in Equation (7), the larger $\alpha$ results in a higher chance of developing the new continuous cluster and adding to the collection of the existing discrete clusters. The forms of the prior and posterior densities used to simulate the main parameters $({\mathit{\theta }}_j^{*},{\alpha }^{*},{\mathit{w}}_j^{*})$ on lines from 16 to 23 in Algorithm (A2) are presented in Appendix A.

There is a couple of points to note. The Gibbs sampler for the DPM described here can be characterized by the use of infinite clusters and covariates. Due to the infinite mixture capacity, the resulting clusters can be kept as homogeneous as possible. In this process, the within-class heterogeneity can be captured between parameters across the observations, and the DPM utilizes such dependencies within existing clusters to determine the rationale for the development of new clusters. The DPM harnesses the power of the covariate as well. For example, the DPM associates individual policies with the unobserved claim (in new clusters) and the observed claims (in old clusters), matching on the covariate information. The investigation of the infinite clusters, covariates, and the continuous cluster selection process in the DPM are briefly illustrated in the diagram in Figure 2. As a result, the unobserved claim problem mentioned in Figure 1 can be addressed by the new cluster introduction, which leads to a better approximation of $S_h$.

3.4. Modelling $S_h|X_h$ with MAR Covariate

The DPM model for complete case data ($S_h,{\mathit{X}}_h)$ has been discussed in Section 3.3. In this Section, we present our novel imputation strategy for the MAR covariate in the DPM framework in which the missing values are explained by the observed data. We focus on the missingness in the binary type covariate. In addition, we specify here different prior distributions and the corresponding posterior distributions constructed for the Gibbs sampler, taking into account the MAR covariate. With the model definition in Equation (1), suppose the binary covariate ${\mathit{x}}_1$ has missingness within it. To handle this MAR covariate, we consider the following modifications in the DPM Gibbs sampler:

a)

Imputation: The missing covariate impacts on the parameter - $\mathit{\theta },\mathit{w}$ - update. For ${\mathit{w}}_j$, only the observations $S_h$ without the missing covariate are used to update. If the cluster does not have any observations with complete data for that covariate, then a draw from the prior distribution would be used to update. For ${\mathit{\theta }}_j$, however, we must first impute values for the missing covariates $x_{1h}$ for all observations $S_h$ within the cluster. Since having already defined a full joint model - $f\left(S_h\right|{\mathit{X}}_h,{\mathit{\theta }}_j)·f\left({\mathit{X}}_h\right|{\mathit{w}}_j)$ - in Section 3.2, we can obtain draws for the MAR covariate $x_{1h}$ from the imputation model such as $f_{Bern}\left(x_{1h}\right|S_h,{\mathit{\theta }}_j,{\mathit{w}}_j)\propto f\left(S_h\right|{\mathit{X}}_h,{\mathit{\beta }}_j,{\sigma }_j^2,{\xi }_j)·f_{Bern}\left(x_{1h}\right|{\pi }_j)$ at each iteration of the Gibbs sampler. The imputation process is briefly illustrated in Figure 3. Once all missing data in all covariates has been imputed, then we can sample from the posterior for $\mathit{\theta }$ and the parameters of each cluster ${\mathit{\beta }}_j,{\sigma }_j^2$ are re-calculated. After this cycle is complete, the imputed data is discarded and the same imputation steps are repeated every iteration.

b)

Re-clustering: To determine each cluster probability after the imputations, the algorithm re-defines the two main components for the cluster probability calculation - 1) covariate model, 2) outcome model. For the covariate model $f\left({\mathit{X}}_h\right|{\mathit{w}}_j)$, we set this equal to the density functions of only those covariates with complete data for observation h. Assuming that ${\mathit{X}}_h=\{x_{1h},x_{2h}\}$, and the covariate ${\mathit{x}}_1$ is missing for observation h, then we drop $x_{1h}$ and only use $x_{2h}$ in the covariate model,

$$f\left({\mathit{X}}_h\right|{\mathit{w}}_j)=f_N\left(x_{2h}\right|{\mathit{w}}_{2j})$$

(8)

This is the refined covariate model for the cluster j with the observation h where the data in ${\mathit{x}}_1$ is not available. For the outcome model $f\left(S_h\right|{\mathit{X}}_h,{\mathit{\theta }}_j)$, the algorithm simply takes the imputation model for each cluster and integrates them out the covariates with missing data. This reduces the degree of variances introduced by the imputations. In our case, as covariate ${\mathit{x}}_1$ is missing for observation h, this missing covariate can be removed from the ${\mathit{X}}_h$ term that is being conditioned on. Therefore, the refined outcome model is

$$f\left(S_h\right|x_{2h},{\mathit{\theta }}_j)\propto \int f\left(S_h\right|{\mathit{X}}_h,{\mathit{\theta }}_j)·f_{Bern}\left(x_{1h}\right|{\mathit{w}}_{1j})d{x}_{1h}$$

(9)

A similar process is conducted for each observation with missing data and each combination of missing covariates. Hence, using Equation (8,9), the cluster probabilities and the predictive distribution can be obtained as illustrated in Step III in Figure 4.

c)

Parameter update: The cluster probability computation is followed by the parameter re-estimation for each cluster, which is illustrated via the diagram in Figure 5. This is the same idea as what we have discussed about the parameter - $\mathit{\theta },\mathit{w}$ - update in Figure 2.

4. Bayesian Inference for $S_h|X_h$ with MAR Covariate

The efficient simulation for the model parameters - $\mathit{\theta }:\{\mathit{\beta },{\sigma }^2,\xi \tilde{\mathit{\beta }}\},\mathit{w}:\{\pi ,\mu ,{\tau }^2\}$, and $\alpha$ - requires the proper parameterization in the parameter models - prior parameter model and posterior parameter model. The accurate estimations of cluster probabilities rely on the legitimate development of data models - outcome model and covariate model - and the model parameter simulation results that govern the data model behaviors. This section is centered on the novel development of parameter and data models, providing the details of the DPM implementation integrated with the MAR imputation strategy.

4.1. Parameter Models and MAR Covariate

Our study is based on a three-level hierarchical structure: the first level regards the data models such as the log skew-normal outcome model and the Bernoulli, Gaussian covariate models, the second level involves the parameter models such as $p\left(\mathit{\theta }\right|S_h,{\mathit{X}}_h),p\left(\mathit{w}\right|{\mathit{X}}_h)$ to explain the data, and the third level is developed from the generalized regression to explain the parameters or the related hyperparameters such as $a_0,b_0,{\nu }_0,c_0,d_0,{\mu }_0,{\tau }_0^2,e_0,{\gamma }_0,g_0$ and $h_0$ to set a probabilistic distribution on the parameter vectors $\mathit{\theta }=\{\mathit{\beta },{\sigma }^2,\xi ,\tilde{\mathit{\beta }}\},\mathit{w}=\{\pi ,\mu ,{\tau }^2\}$. See Variable Definition for further information on the variables. Given the model definition in Equation (1), we consider a set of conjugate parameter models due to its computational advantages [21]. For $S_h\sim \delta \left({\mathit{X}}_h^T{\tilde{\mathit{\beta }}}_j\right)\mathbb{1}(S_h=0)+\left[1-\delta \left({\mathit{X}}_h^T{\tilde{\mathit{\beta }}}_j\right)\right]LogSN({\mathit{X}}_h^T{\mathit{\beta }}_j,{\sigma }_j^2,{\xi }_j)$, ${\mathit{x}}_1\sim Bern\left({\pi }_j\right)$, and ${\mathit{x}}_2\sim N({\mu }_j,{\tau }_j^2)$, the prior models come in

$$\begin{array}{cccccc}\hfill p_0\left({\sigma }_j^2\right|a_0,b_0):& \mathit{InvGa}(a_0,b_0),\hfill & \hfill p_0\left({\mathit{\beta }}_j\right|{\mathit{\beta }}_0,{\Sigma }_0):& \mathit{MVN}({\mathit{\beta }}_0,{\sigma }_j^2{\Sigma }_0),\hfill & \hfill & p_0\left({\xi }_j\right|{\nu }_0):\mathit{T}\left({\nu }_0\right)\hfill \\ \hfill p_0\left({\tilde{\mathit{\beta }}}_j\right|{\tilde{\mathit{\beta }}}_0,{\tilde{\Sigma }}_0):& \mathit{MVN}({\tilde{\mathit{\beta }}}_0,{\tilde{\Sigma }}_0),\hfill & \hfill p_0\left({\pi }_j\right|c_0,d_0):& \mathit{Beta}(c_0,d_0),\hfill & \hfill & p_0\left({\mu }_j\right|{\mu }_0,{\tau }_0^2):\mathit{N}({\mu }_0,{\tau }_j^2),\hfill \\ \hfill p_0\left({\tau }_j^2\right|e_0,{\gamma }_0):& \mathit{InvGa}(e_0,{\gamma }_0),\hfill & \hfill p_0\left(\alpha \right|g_0,h_0):& \mathit{Ga}(g_0,h_0)\hfill \end{array}$$

and their corresponding kernels chosen in this study are listed in Appendix A.1. Accordingly, the Dirichlet process prior (probability measure) $G_0$ in our case can be defined as $G_0=\mathit{MVN}({\mathit{\beta }}_0,{\Sigma }_0)\times \mathit{InvGa}(a_0,b_0)\times \mathit{T}\left({\nu }_0\right)\times \mathit{MVN}({\tilde{\mathit{\beta }}}_0,{\tilde{\Sigma }}_0)\times \mathit{Beta}(c_0,d_0)\times \mathit{N}({\mu }_0,{\tau }_j^2)\times \mathit{InvGa}(e_0,{\gamma }_0)\times \mathit{Ga}(g_0,h_0)$. With a feed of the observed data inputs - $(S_h,x_{1h},x_{2h})$ -, the prior models for each cluster j described above will be updated into the following posterior models analytically apart from ${\mathit{\theta }}_j=\{{\mathit{\beta }}_j,{\sigma }_j^2,{\xi }_j,{\tilde{\mathit{\beta }}}_j\}$.

$$\begin{array}{cc}\hfill & p\left({\pi }_j\right|c_0,d_0,S,{\mathit{x}}_1):\mathit{Beta}(c_{new},d_{new})\hfill \\ \hfill & p\left({\mu }_j\right|{\mu }_0,{\tau }_0^2,S,{\mathit{x}}_2):\mathit{N}({\mu }_{new},{\tau }_{new}^2),p\left({\tau }_j^2\right|e_0,{\gamma }_0,S,{\mathit{x}}_2):\mathit{InvGa}(e_{new},{\gamma }_{new})\hfill \\ \hfill & p\left(\alpha \right|g_0,h_0,h,J,\eta ,{\pi }_{\eta }):{\pi }_{\eta }\mathit{Ga}(g_0+J,h_0-log\left(\eta \right))+(1-{\pi }_{\eta })\mathit{Ga}(g_0+J-1,h_0-log\left(\eta \right))\hfill \end{array}$$

(10)

and their corresponding parameterizations are elaborated in Appendix A.2. Note that the value of the precision parameter $\alpha$ relies on the total cluster number J, thus does not vary by the cluster membership j, and its derivation of the posterior parameterization is not subject to the Bayesian conjugacy. Hence, we instead adapt the form of the posterior density for the $\alpha$ suggested by Escobar and West (1995) [22], and its derivation is shown in Appendix C.1. As for ${\mathit{\theta }}_j=\{{\mathit{\beta }}_j,{\sigma }_j^2,{\xi }_j,{\tilde{\mathit{\beta }}}_j\}$, there are no conjugate priors available for log skew-normal likelihood, but their posterior samples can be secured by the conventional metropolis hastings described in Algorithm (A2) in Appendix A.

Considering that ${\mathit{x}}_1$ has missing data, although the parameterizations of the posterior densities for the covariate parameter model of $\mathit{w}$ and the precision $\alpha$ listed in Equation (10) are not affected, any outcome data of $S_h$ with missingness should be dropped; therefore, $n_j$ and ${\mathit{x}}_1$ are defined with the only observations in cluster j that are not missing. This imputation example is provided in Appendix C.2. For the outcome parameter model of ${\mathit{\theta }}_j$, the missing covariate ${\mathit{x}}_1$ must be imputed before its posterior computation shown in Algorithm (A2). Once the parameters are updated with the imputation, the data models can be constructed as described in Equation (8,9).

4.2. Data Models and MAR Covariate

Data models are the main components for cluster probability computations depicted in Figure 2. As with the development of parameter models, the covariate data model of $\mathit{X}$ ignores the observations with missingness while the outcome data model of $S_h$ requires to complete the covariates beforehand. However, the formulation of their densities can be more complex due to the marginalization process with respect to the missing covariate. In addition, as discussed in Section 3.2, the data model development is bound by the types of clusters such as discrete clusters $f\left(S_h\right|{\mathit{X}}_h,{\mathit{\theta }}_j),f\left({\mathit{X}}_h\right|{\mathit{w}}_j)$ and continuous clusters $f_0\left(S_h\right|{\mathit{X}}_h),f_0\left({\mathit{X}}_h\right)$.

a)

covariate model for the discrete cluster: $f\left({\mathit{X}}_h\right|{\mathit{w}}_j)$

Focusing on the scenario that ${\mathit{x}}_1$ is binary, ${\mathit{x}}_2$ is Gaussian, and the only covariate with missingness is $x_{1h}$, we simply drop the covariate $x_{1h}$ to develop the covariate model for the discrete cluster. For instance, when computing the covariate probability term for hth observation in j cluster, the covariate model $f(x_{1h},x_{2h}|{\pi }_j,{\mu }_j,{\tau }_j^2)$ simply becomes $f\left(x_{2h}\right|{\mu }_j,{\tau }_j^2)$ due to the missingness of $x_{1h}$. As we have ${\mathit{x}}_2$ that is assumed to be normally distributed as defined in Equation (1), its probability term is

$$f\left(x_{2h}\right|{\mu }_j,{\tau }_j^2)=\frac{1}{\sqrt{2\pi {\tau }_j^2}}exp\left\{\frac{-{(x_{2h}-{\mu }_j)}^2}{2{\tau }_j^2}\right\}$$

(11)

instead of

$$f(x_{1h},x_{2h}|{\pi }_j,{\mu }_j,{\tau }_j^2)={\pi }_j^{x_{1h}}{\left(1-{\pi }_j\right)}^{1-x_{1h}}·\frac{1}{\sqrt{2\pi {\tau }_j^2}}exp\left\{\frac{-{(x_{2h}-{\mu }_j)}^2}{2{\tau }_j^2}\right\}$$

b)

covariate model for the continuous cluster: $f_0\left({\mathit{X}}_h\right)$

If the binary covariate $x_{1h}$ is missing, by the same logic, we drop the covariate $x_{1h}$ for the continuous cluster; however, using Equation (4), the covariate model for the continuous cluster integrates out the relevant parameters simulated from the Dirichlet process prior $G_0$ as follows:

$$\begin{array}{cc}\hfill f_0\left(x_{2h}\right)=\int f(x_{2h}& |\mu ,{\tau }^2)d{G}_0(\mu ,{\tau }^2)=\int f\left(x_{2h}\right|\mu ,{\tau }^2)·p\left(\mu \right|{\tau }^2)·p\left({\tau }^2\right)d\mu d{\tau }^2\hfill \\ \hfill & =\frac{{\gamma }_0^{e_0}\Gamma (e_0+1/2)}{2\sqrt{\pi }\Gamma \left(e_0\right)}{\left({\gamma }_0+\frac{{(x_{2h}-{\mu }_0)}^2}{4}\right)}^{-(e_0+1/2)}\hfill \end{array}$$

(12)

instead of

$$\begin{array}{cc}\hfill f_0(& x_{1h},x_{2h})=\int f(x_{1h},x_{2h}|\pi ,\mu ,{\tau }^2)·p\left(\pi \right)·p\left(\mu \right|{\tau }^2)·p\left({\tau }^2\right)d\pi d\mu d{\tau }^2\hfill \\ \hfill & =\frac{\mathit{B}(x_{1h}+c_0,1-x_{1h}+d_0)}{\mathit{B}(c_0,d_0)}·\frac{{\gamma }_0^{e_0}\Gamma (e_0+1/2)}{2\sqrt{\pi }\Gamma \left(e_0\right)}{\left({\gamma }_0+\frac{{(x_{2h}-{\mu }_0)}^2}{4}\right)}^{-(e_0+1/2)}\hfill \end{array}$$

The derivation of the distributions above is provided in Appendix C.3.

c)

outcome model for the discrete cluster: $f\left(S_h\right|{\mathit{X}}_h,{\mathit{\theta }}_j)$

In developing the outcome model, as with the parameter model case discussed in Section 4.1 and Appendix C.2, it should be ensured that the covariate is complete beforehand. With all missing data in $x_{1h}$ imputed, the outcome model for the discrete cluster is obtained by marginalizing the joint - $f(S_h,x_{1h}|x_{2h},{\mathit{\theta }}_j,{\pi }_j)$ - out the MAR covariate $x_{1h}$, which is a log skew-normal mixture as follows:

$$\begin{array}{cc}& f\left(S_h\right|x_{2h},{\mathit{\beta }}_j,{\sigma }_j^2,{\xi }_j,{\tilde{\mathit{\beta }}}_j)=\sum _{x_{1h}=0}^{1}f\left(S_h\right|x_{1h},x_{2h},{\mathit{\beta }}_j,{\sigma }_j^2,{\xi }_j,{\tilde{\mathit{\beta }}}_j)·f\left(x_{1h}\right|{\pi }_j)\hfill \\ \hfill & =f(S_h,x_{1h}=1|x_{2h},{\mathit{\beta }}_j,{\sigma }_j^2,{\xi }_j,{\tilde{\mathit{\beta }}}_j,{\pi }_j)+f(S_h,x_{1h}=0|x_{2h},{\mathit{\beta }}_j,{\sigma }_j^2,{\xi }_j,{\tilde{\mathit{\beta }}}_j,{\pi }_j)\hfill \\ \hfill & =\delta \left({\mathit{X}}_h^T{\tilde{\mathit{\beta }}}_j\right)\mathbb{1}(S_h=0)+\left[1-\delta \left({\mathit{X}}_h^T{\tilde{\mathit{\beta }}}_j\right)\right]·{\displaystyle \frac{2}{{\sigma }_j{S}_h}}\hfill \\ \hfill & ·\varphi \left(\frac{log{S}_h-({\beta }_{j0}+{\beta }_{j1}+{\beta }_{j2}{x}_{2h})}{{\sigma }_j}\right)·\mathrm{\Phi }\left({\xi }_j\frac{log{S}_h-({\beta }_{j0}+{\beta }_{j1}+{\beta }_{j2}{x}_{2h})}{{\sigma }_j}\right){\pi }_j\hfill \\ \hfill & +\delta \left({\mathit{X}}_h^T{\tilde{\mathit{\beta }}}_j\right)\mathbb{1}(S_h=0)+\left[1-\delta \left({\mathit{X}}_h^T{\tilde{\mathit{\beta }}}_j\right)\right]·{\displaystyle \frac{2}{{\sigma }_j{S}_h}}\hfill \\ \hfill & ·\varphi \left(\frac{log{S}_h-({\beta }_{j0}+{\beta }_{j2}{x}_{2h})}{{\sigma }_j}\right)·\mathrm{\Phi }\left({\xi }_j\frac{log{S}_h-({\beta }_{j0}+{\beta }_{j2}{x}_{2h})}{{\sigma }_j}\right)·(1-{\pi }_j)\hfill \end{array}$$

(13)

instead of

$$\begin{array}{cc}\hfill & f\left(S_h\right|x_{1h},x_{2h},{\mathit{\beta }}_j,{\sigma }_j^2,{\xi }_j,{\tilde{\mathit{\beta }}}_j)\hfill \\ \hfill & =\delta \left({\mathit{X}}_h^T{\tilde{\mathit{\beta }}}_j\right)\mathbb{1}(S_h=0)+\left[1-\delta \left({\mathit{X}}_h^T{\tilde{\mathit{\beta }}}_j\right)\right]·{\displaystyle \frac{2}{{\sigma }_j{S}_h}}\hfill \\ \hfill & ·\varphi \left(\frac{log{S}_h-({\beta }_{j0}+{\beta }_{j1}{x}_{1h}+{\beta }_{j2}{x}_{2h})}{{\sigma }_j}\right)·\mathrm{\Phi }\left({\xi }_j\frac{log{S}_h-({\beta }_{j0}+{\beta }_{j1}{x}_{1h}+{\beta }_{j2}{x}_{2h})}{{\sigma }_j}\right)\hfill \end{array}$$

d)

outcome model for the continuous cluster: $f_0\left(S_h\right|{\mathit{X}}_h)$

Once a missing covariate ${\mathit{x}}_1$ is fully imputed and the outcome model is marginalized out conditioned to the MAR covariate $x_{1h}$, the outcome model $f_0\left(S_h\right|x_{2h})$ for the continuous cluster can also be computed by integrating out the relevant parameters, using Equation (4).

$$f_0\left(S_h\right|x_{2h})=\int f\left(S_h\right|x_{2h},\mathit{\beta },{\sigma }^2,\xi ,\tilde{\mathit{\beta }})·p\left(\mathit{\beta }\right)·p\left({\sigma }^2\right)·p\left(\xi \right)·p\left(\tilde{\mathit{\beta }}\right)d\mathit{\beta }d{\sigma }^{2}d\xi d\tilde{\mathit{\beta }}$$

(14)

However, it can be too complicated to compute its form analytically. Instead, we can integrate the joint model out the parameters, using Monte Carlo integration. For example, we can do the following for each $h=1,\cdots ,H$.

(i)

Sample $\mathit{\beta },{\sigma }^2,\xi ,\tilde{\mathit{\beta }}$ from the DP prior densities $G_0$ specified previously.

(ii)

Plug in these samples into $f\left(S_h\right|x_{2h},\mathit{\beta },{\sigma }^2,\xi ,\tilde{\mathit{\beta }})·p\left(\mathit{\beta }\right)·p\left({\sigma }^2\right)·p\left(\xi \right)·p\left(\tilde{\mathit{\beta }}\right)$.

(iii)

Repeat the above steps many times, recording each output.

(iv)

Divide the sum of all output values by the number of Monte Carlo samples, which will be the approximate integral.

4.3. Gibbs Sampler Modification for MAR Covariate

We have examined the parameter models and data models to update the parameters of the DPM based on probabilistically imputed values of the MAR covariate. Now we set out some modifications of the DPM and let the Gibbs sampler in Algorithm (A2) in Appendix B. address the MAR covariate of ${\mathit{x}}_1$. The Gibbs sampler will alternate between imputing missing data and drawing parameters until it reaches a stationary distribution of the parameters. We elaborate below on the modifications that fit into Algorithm (A2) to update the clustering scenarios and the posterior cluster parameters properly.

a)

In line 6, with the presence of missing covariate $x_{1h}$, the modification of the cluster probability for the observation $(S_h,x_{1h},x_{2h})$ that belongs to the discrete cluster j can be made as follows,

$$P(s_h=j)=p\left(s_h\right|s_{-h})·f\left(x_{2h}\right|{\mu }_j,{\tau }_j^2)·f\left(S_h\right|x_{2h},{\mathit{\beta }}_j,{\sigma }_j^2,{\xi }_j,{\tilde{\mathit{\beta }}}_j)$$

where $f\left(x_{2h}\right|{\mu }_j,{\tau }_j^2)$ is from Equation (11), and $f\left(S_h\right|x_{2h},{\mathit{\beta }}_j,{\sigma }_j^2,{\xi }_j,{\tilde{\mathit{\beta }}}_j)$ is from Equation (13).

b)

In line 9, with the presence of missing covariate $x_{1h}$, the modification of the cluster probability for the observation $(S_h,x_{1h},x_{2h})$ that belongs to the continuous cluster $J+1$ can be made as follows,

$$P(s_h=J+1)=p\left(s_h\right|s_{-h})·f_0\left(x_{2h}\right)·f_0\left(S_h\right|x_{2h})$$

where $f_0\left(x_{2h}\right)$ is from Equation (12), and $f_0\left(S_h\right|x_{2h})$ is from Equation (14).

c)

In line 22, with the presence of missing covariate $x_{1h}$, the imputation should be made before simulating the parameter ${\mathit{\theta }}_j^{*}$ as follows,

The imputation model formulation in the above has been discussed in Section 3.4.

Again, these modifications allow to draw missing covariate values from the conditional posterior density at each iteration, using the Metropolis-Hastings with a random walk.

5. Empirical Study

5.1. Data

The performance of our DPM framework is assessed based on two insurance datasets. They highlight data difficulties such as unobservable heterogeneity in an outcome variable and MAR covariates. For simplicity, in each dataset, we only consider two covariates - one binary and one continuous - to explain its loss information (outcome variable). In this study, all computations on these two datasets are performed in the same data format:

The first dataset is PnCdemand, which is about the international property and liability insurance demand of 22 countries over 7 years from 1987 to 1993. Secondly, we use a dataset drawn from the Wisconsin Local Government Property Insurance Fund (LGPIF) with information about the insurance coverage for government building units in Wisconsin for years from 2006 to 2010. The first one - PnCdemand - can be obtained from the R package CASdatasets. The dataset is relatively small as it has $H=240$ cases with an outcome variable GenLiab: the individual loss amount under the policies of general insurance for each case. As for covariates, we consider one indicator variable of the statutory law system (LegalSyst:1 or 0) and one continuous variable that measures a risk aversion rate (RiskAversion) for each area. For additional background on this dataset, see Browne et al. (2000)  [23]. In the LGPIF dataset, the insurance coverage samples for the government properties from $H=5660$ policies are provided. The outcome variable is the sum of all types of losses (Total Losses) for each policy. Only the covariates - LnCoverage, Fire5 - are considered in our study. Fire5 is a binary covariate that indicates fire-protection levels while LnCoverage is a continuous covariate that informs a total coverage amount in a logarithmic scale. For further details, see Quan et al. [24].

Histograms of the losses of the two datasets are exhibited in Figure 6. Due to the significant skewness, the loss data are log-transformed to attain Gaussianity. As shown in the histograms, each distribution displays different characteristics in regard to skewness, modality, excess of zeros, etc. Note that the zero-inflated outcome variable in LGPIF data (b1, b2 in Figure 6) requires a two-part modeling technique that distinguishes the probabilities of the outcome being zero and positive.

5.2. Three Competitor Models and Evaluation

Our DPM framework is compared to other commonly used actuarial models in practice. We employ three predictive models as benchmarks - namely, a generalized linear mixture model (GLM), multivariate adaptive regression spline (MARS), and generalized additive model (GAM). In each dataset, we assume different distributions for the outcome variables, and thus the three benchmark models are built upon the different outcome data models. For example, the PnCdemand dataset (a1,a2) that appeared in Figure 6, has a high frequency of small losses without zero values, hence it is safe to use a gamma mixture to explain the outcome data. As for the LGPIF data (b1,b2) in Figure 6, we consider the outcome data model based on a Tweedie distribution to accommodate the zero-inflated loss data. The benchmark models are implemented in R with the mgcv, splines, and mice packages.

All four models are trained, and investigations are performed in terms of model fit, prediction accuracy, and the conditional tail expectation (CTE) of the predictive distribution. Note that the goodness of fit value for a DPM is not available in Table 1,Table 2. Teh (2010) [25] argues that the goodness of fit evaluation for a DPM is unnecessary as underfitting is mitigated by the unbounded complexity of a DPM while overfitting is alleviated by the approximation of posterior densities over each parameter in a DPM. Gelman et al. (2007) [26] point out Posterior predictive check, which compares the simulated data under the fitted DPM to the observed data, can be useful in studying model adequacy, but its usage cannot be for model comparison. Therefore, the goodness of fit is only compared between the rival models. For the evaluation of prediction performance, the sum of square prediction error (SSPE) and sum of square absolute error (SAPE) are used.

5.3. Result 01. International General Insurance Liability Data

For this dataset, a training set of response and covariates pair $(Y,\mathit{X})$ with n = 160 records, and a test set of response and covariates pair $(Y^{\prime },{\mathit{X}}^{\prime })$ with m = 80 records are constructed. We implement the following DPM:

$$\begin{array}{cc}\hfill Y_h|{\mathit{x}}_1,{\mathit{x}}_2,{\mathit{\beta }}_j,{\sigma }_j^2& \sim LogN\left({\mathit{X}}^T{\mathit{\beta }}_j,{\sigma }_j^2\right)\hfill \\ \hfill {\mathit{x}}_1|{\pi }_j& \sim Bern\left({\pi }_j\right)\hfill \\ \hfill {\mathit{x}}_2|{\mu }_j,{\tau }_j^2& \sim N\left({\mu }_j,{\tau }_j^2\right)\hfill \\ \hfill \{{\mathit{\theta }}_j,{\mathit{w}}_j\}& \sim G\hfill \\ \hfill G& \sim DP\left(\alpha ,G_0\right)\hfill \end{array}$$

A log-normal likelihood is chosen to accommodate the individual loss $Y_h$:GenLiab for a policy h. The covariate $x_2$:RiskAversion is subject to missingness, and found to depend on $Y_h$ (a MAR case). This is addressed by the internalized imputation process as discussed in Figure 3. The posterior parameters of ${\mathit{\theta }}_j,{\mathit{w}}_j$ are estimated with our DPM Gibbs sampler presented in Algorithm (A2). The algorithm runs 10,000 iterations until convergence, and the resulting scenarios of clustering mixture are shown in Figure 7. The plot reveals the overlays of predictive densities on the log scale from the last 100 iterations that are tied to convergence.

Figure 8 lists the classical data imputation process - Multivariate Imputation Chained Equation (MICE) - and predictive densities produced from our rival models - GLM, GAM, MARS. The MICE runs multiple imputation chains, and selects the imputation values from the final iteration. This process results in multiple candidate datasets. The trace plots (a1,a2) monitor the imputation mean and variance for the missing values in the dataset. In the covariate distribution plot (a3), the density of the observed covariate shown in blue is compared with the ones of the imputed covariate for each imputed dataset shown in red. The parameter inferences for the rival models are performed based on the imputed datasets tied to convergence [27]. The gamma distribution is chosen to fit the rival models as the $Y_h$ is continuous and positively skewed with a constant coefficient of variation. The gamma-based predictive density plots (b1,b2,b3) estimated with GLM, GAM, MARS look similar, showing unusual bumps near the right tail.

In Figure 9, a histogram of the outcome data in the test set is displayed. The posterior mean densities for out-of-sample predictions produced with our DPM along with the rival models’ density estimates are overlaid on the histogram. Judging from the plot, one can say that our DPM model generates the best approximation. While the rival models generate smooth, mounded curves to make predictions, our DPM captures all possible peaks and bumps, which is closer to the actual situation. According to Table 1, the rival models produce slightly higher SAPEs, but lower SSPEs, compared to our proposed DPM. As SAPE weights all the individual differences equally, we can assume that the rival models tend to give too much focus on the most probable data points and miss some outliers. This is mainly due to the insufficient sample size. However, our DPM has good performance under small sample sizes when there is sufficient prior knowledge available. From the perspective of CTE, Table 1 shows that our DPM proposes a heavier tail than other rival models, which reflects that our DPM captures more uncertainties given the small sample size.

5.4. Result 02. LGPIF Data

For this dataset, a training set of response and covariates pair $(S,\mathit{X})$ with n = 4529 records, and a test set of response and covariates pair $(S^{\prime },{\mathit{X}}^{\prime })$ with m = 1110 records are constructed. We implement the following DPM:

$$\begin{array}{cc}\hfill S_h|{\mathit{x}}_1,& {\mathit{x}}_2,{\mathit{\beta }}_j,{\sigma }_j^2,{\xi }_j,{\tilde{\mathit{\beta }}}_j\hfill \\ \hfill & \sim \delta \left({\mathit{X}}^T{\tilde{\mathit{\beta }}}_j\right)\mathbb{1}(S_h=0)+\left[1-\delta \left({\mathit{X}}^T{\tilde{\mathit{\beta }}}_j\right)\right]LogSN\left({\mathit{X}}^T{\mathit{\beta }}_j,{\sigma }_j^2,{\xi }_j\right)\hfill \\ \hfill {\mathit{x}}_1|{\pi }_j& \sim Bern\left({\pi }_j\right)\hfill \\ \hfill {\mathit{x}}_2|{\mu }_j,{\tau }_j^2& \sim N\left({\mu }_j,{\tau }_j^2\right)\hfill \\ \hfill \{{\mathit{\theta }}_j,{\mathit{w}}_j\}& \sim G\hfill \\ \hfill G& \sim DP\left(\alpha ,G_0\right)\hfill \end{array}$$

As the outcome $S_h$:Total Losses for a policy h in this dataset is considered to be distributed with the sum of log-normal densities, a log skew-normal likelihood is chosen to approximate this convolution [7]. The covariate $x_1$:Fire5 is subject to missingness under MAR, and the internalized imputation process illustrated in Figure 3 resolves this issue without creating imputed datasets. As the outcome $S_h$ exhibits zero inflation, we employ a two-part model, using a sigmoid and indicator function. Our DPM Gibbs sampler described in Algorithm (A2) produces the posterior parameters of ${\mathit{\theta }}_j,{\mathit{w}}_j$ with 10,000 iterations until convergence. Figure 10 reveals the resulting scenarios of clustering mixture. In the plot, there are 100 predictive densities suggested by our DPM, each of which stands for the convergence of the estimation results.

The output of the MICE and the resulting predictive densities from the rival models are displayed in Figure 11. The rival models are built upon a Tweedie distribution due to its ability to account for a large number of zero losses, and the flexibility to capture the unique loss patterns of the different classes of policyholders. According to the plot, all three rival models reasonably capture zero inflation, but the GAM tends to suggest more bumps that indicate a need for further assessment of the prediction uncertainty.

The overall out-of-sample prediction comparison is made in the histogram overlayed with predictive density curves generated from the four models in Figure 12. From the plot, it is apparent that the posterior predictive density proposed by our DPM best explains the new samples while other rival models keep producing multiple peaks. The prediction performance of our DPM is confirmed by the smallest SSPE and SAPE in Table 2. In terms of CTE, all three rival models suggest a similar level of tailedness, reflecting the knowledge obtained from the observed data. However, our DPM goes beyond this and proposes a much heavier tail. This is because our DPM accommodates the presence of outliers and shapes the tail behavior based on the combined knowledge of prior parameters and the observations available.

6. Discussion

This paper proposes a novel DPM framework for actuarial practice to model total losses with the incorporation of MAR covariates. Both log-normal and log skew-normal DPM present overall good empirical performances in capturing the shape of the distribution, out-of-sample prediction, and the estimation of the tailedness. This suggests that it is worth considering our DPM framework in order to avoid various model risks or biases in insurance claim analysis.

6.1. Research Questions

Regarding RQ1, we propose a DPM framework to address the within-cluster heterogeneity emerging from the inclusion of covariates. By allowing for an infinite number of clustering scenarios determined by the observations as well as prior knowledge, our DPM outperforms the rival methods in drawing the lines for the cluster membership. This can be assessed by examining the homogeneity of the resulting clusters. In our case, we fit cluster-wise GLMs (based on Gamma and Tweedie) to the data points within each resulting cluster to compare the goodness-of-fit, and the consistent AICs across all clusters endorse the benefits of the DPM. Similarly, our rival methods such as GAM or MARS can capture heterogeneity by using customized smooth functions across different subsets of the data, but we observe some statistically insignificant smooth terms, indicating the presence of heterogeneity in the cluster.

In terms of RQ2, we suggest incorporating the imputation steps into the parameter and cluster membership update process in the DPM Gibbs sampler by leveraging the joint distribution of the observed outcomes and missing covariates. This approach allows the imputed values to be consistent with the observed data, preserving the correlation structure within the dataset. In order to make a comparison of our approach with an existing alternative, we additionally employ a chained equation technique. The multiple sets of imputed values simulated from both approaches are investigated, and the result shows that our DPM Gibbs sampler does not represent a significant improvement over the chained equation because their average estimates of the imputed values are closer to each other. However, we feel that this result is mainly due to the relatively low dimensionality of the datasets we use and their simple data structure. The specific characteristics or dependencies in the data and the complexity of the missing patterns would give different results in practice.

As for RQ3, we fit a log skew-normal density to the aggregate loss outcomes. In order to assess its performance, one can consider Minimax approximation, Least squares approximation, Log shifted gamma approximation, etc. as the competitors. Li (2008) [7] provides a useful comparison between these competitors by overlaying the cumulative density curves for each technique, but its experiments are grounded on the simulated log-normal data with the pre-defined parameters and assumptions, which cannot be easily controlled in real-world scenarios. Therefore, we feel that the choice of the best approximation technique should be made based on the identification of the specific characteristics of the dataset. In our case, each summand in our dataset is significantly different from each other in magnitudes (the Minimax is inappropriate) and LGPIF data has a large volume of data smaller than 5 (the Log shifted gamma is inappropriate); therefore, we choose a log skew-normal density that is relatively simple while giving an accurate approximation at the lower region of the distribution.

6.2. Future Work

There are several concerns with our log skew-normal DPM framework.

(a)

Dimensionality: First, in our analysis, we only use two covariates (binary and continuous) for simplicity, hence more complex data should be considered. As the number of covariates grows, the likelihood components (covariate models) to describe the covariates grow, which results in the shrinking of the cluster weights. Therefore, using more covariates might enhance the level of sensitivity and accuracy in the creation of cluster memberships. However, it can also introduce more noise or hidden structures that render the resulting predictive distributions unstable. In this sense, further research on the problem of high dimensional covariates in the DPM framework would be worthwhile.

(b)

Measurement error: Second, although our focus in this article is MAR covariate, mismeasured covariate is an equally significant challenge that impairs the proper model development in insurance practice. For example, Aggarwal et al. (2016) [28] point out that “model risk" mainly arises due to missingness and measurement error in variables, leading to flawed risk assessments and decision-making. Thus, further investigation is necessary to explore the specialized construction of the DPM Gibbs sampler for mismeasured covariates, aiming to prevent the issue of model risk.

(c)

Sum of log skew-normal: Third, as an extension to the approximation of total losses $S_h$ (the sum of individual losses) for a policy, we recommend researching into ways to approximate the sum of total losses $\tilde{S}$ across entire policies. In other words, we pose the question of “how to approximate the sum of log skew-normal random variables". From the perspective of an executive or an entrepreneur whose concern is the total cash flow of the firm, nothing might be more important than the accurate estimation of the sum of total losses in order to identify the insolvency risk or to make important business decisions.

(d)

Scalability: Lastly, we suggest investigating the scalability of the posterior simulation by our DPM Gibbs sampler. As shown in our empirical study on the PnCdemand dataset, our DPM framework produces reliable estimates with relatively small sample sizes ($n\le 160$). This is because our DPM framework actively utilizes significant prior knowledge in posterior inference rather than heavily relying on the actual features of the data. In the result from the LGPIF dataset, our DPM exhibits stable performance at sample size $n=4529$ as well. However, a sample size of over 10000 is not explored in this paper. With increasing amounts of data, our DPM framework raises a question of computational efficiency due to the growing demand for computational resources or degradation in performance [29]. This is an important consideration, especially in scenarios where the insurance loss information is expected to grow over time.

7. Variable Definitions

    The following variables and functions are used in this manuscript:

$i=1,\dots ,N_h$	observation index i in policy h.
$h=1,\dots ,H$	policy index h with sample (policy) size H.
$j=1,\dots ,J$	cluster index for J clusters.
$s_h$	cluster index $j=1,\dots ,J$ for observation h.
$n_j$	number of observations in cluster j.
$n_j^{-h}$	number of observations in cluster j where observation h removed from.
$Y_{ih}$	individual loss i in a policy observation h.
$S_h$	outcome variable which is a $\Sigma Y_{ih}$ in a policy observation h.
$\tilde{S}$	outcome variable which is a $\Sigma S_h$ across entire policies
${\mathit{X}}_h$	vector of covariates (including ${\mathit{x}}_1,{\mathit{x}}_2$) for observation h.
${\mathit{x}}_1$	vector of covariate (Fire5).
${\mathit{x}}_2$	vector of covariate (Ln(coverage)).
$x_1$	individual value of covariate (Fire5).
$x_2$	individual value of covariate (Ln(coverage)).
$p_0(·)$	parameter model (for prior).
$p(·)$	parameter model (for posterior).
$f_0(·)$	data model (for continuous cluster).
$f(·)$	data model (for discrete cluster).
$\delta (·)$	logistic sigmoid function - expit(·) - to allow for a positive probability of the zero outcome.
${\mathit{\theta }}_j$	set of parameters - $\mathit{\beta },{\sigma }^2,\xi$ - associated with the $f\left(\Sigma Y\right|\mathit{X})$ for j cluster.
${\mathit{w}}_j$	set of parameters - $\pi ,\mu ,\tau$ - associated with the $f\left(\mathit{X}\right)$ for j cluster.
${\mathit{\omega }}_j$	cluster weights (mixing coefficient) for j cluster.
${\mathit{\beta }}_0,{\Sigma }_0$	vector of initial regression coefficients and variance-covariance matrix, i.e. ${\widehat{\sigma }}^2{\left({\mathit{X}}^T\mathit{X}\right)}^{-1}={\mathit{X}}^T\mathit{X}{(\Sigma Y-\Sigma \widehat{Y})}^T(\Sigma Y-\Sigma \widehat{Y})/(n-p)$ obtained from the baseline multivariate Gamma regression of $\Sigma \widehat{Y}>0$.
${\mathit{\beta }}_j$	regression coefficient vector for a mean outcome estimation.
${\sigma }_j^2$	cluster-wise variation value for the outcome.
${\xi }_j$	skewness parameter for log skew-normal outcome.
${\tilde{\mathit{\beta }}}_0,{\tilde{\Sigma }}_0$	vector of initial regression coefficients and variance-covariance matrix obtained from the baseline multivariate logistic regression of $\Sigma \widehat{Y}=0$.
${\tilde{\mathit{\beta }}}_j$	regression coefficient vector for a logistic function to handle zero outcomes.
${\pi }_j$	proportion parameter for Bernoulli covariate.
${\mu }_j,{\tau }_j$	location and spread parameter for Gaussian covariate.
$\alpha$	precision parameter that controls the variance of the clustering simulation. For instance, a larger $\alpha$ allows to select more clusters.
$G_0$	prior joint distribution for all parameters in the DPM - $\beta ,{\sigma }^2,\xi ,\pi ,\mu ,\tau$, and $\alpha$. It allows all continuous, integrable distributions to be supported while retaining theoretical properties and computational tractability such as asymptotic consistency, efficient posterior estimation, etc.
$a_0,b_0$	hyperparameters for Inverse Gamma density of ${\sigma }_j^2$.
$c_0,d_0$	hyperparameters for Beta density of ${\pi }_j$.
${\nu }_0$	hyperparameters for Student’s t density of ${\xi }_j$.
${\mu }_0,{\tau }_0^2$	hyperparameters for Gaussian density of ${\mu }_j$.
$e_0,{\gamma }_0$	hyperparameters for Inverse Gamma density of ${\tau }_j^2$.
$g_0,h_0$	hyperparameters for Gamma density of $\alpha$.
$\eta$	random probability value for Gamma mixture density of the posterior on $\alpha$.
${\pi }_{\eta }$	mixing coefficient for Gamma mixture density of the posterior on $\alpha$.

Appendix A.1. Prior Kernel for Distributions of Outcome, Covariates, and Precision

$$\begin{array}{cc}\hfill p_0\left({\mathit{\beta }}_j\right|{\mathit{\beta }}_0,{\Sigma }_0)& :\mathit{MVN}{({\mathit{\beta }}_0,{\sigma }_j^2{\Sigma }_0)}^{*}\propto e^{\left\{{({\mathit{\beta }}_j-{\mathit{\beta }}_0)}^T{\Sigma }_0^{-1}({\mathit{\beta }}_j-{\mathit{\beta }}_0)\right\}},p_0\left({\sigma }_j^2\right|a_0,b_0):\mathit{InvGa}(a_0,b_0)\propto {\left({\sigma }_j^2\right)}^{-(a_0+1)}·e^{-b_0/{\sigma }_j^2}\hfill \\ \hfill p_0\left({\xi }_j\right|{\nu }_0)& :\mathit{T}\left({\nu }_0\right)\propto {\left(\frac{{\xi }_j^2}{{\nu }_0}+1\right)}^{-({\nu }_0+1)/2},p_0\left({\tilde{\mathit{\beta }}}_j\right|{\tilde{\mathit{\beta }}}_0,{\tilde{\Sigma }}_0):\mathit{MVN}{({\tilde{\mathit{\beta }}}_0,{\tilde{\Sigma }}_0)}^{*}\propto e^{\left\{{({\tilde{\mathit{\beta }}}_j-{\tilde{\mathit{\beta }}}_0)}^T{\tilde{\Sigma }}_0^{-1}({\tilde{\mathit{\beta }}}_j-{\tilde{\mathit{\beta }}}_0)\right\}}\hfill \\ \hfill p_0\left({\pi }_j\right|c_0,d_0)& :\mathit{Beta}(c_0,d_0)\propto {\pi }_j^{(c_0-1)}·{(1-{\pi }_j)}^{(d_0-1)},p_0\left({\mu }_j\right|{\mu }_0,{\tau }_0^2):\mathit{N}({\mu }_0,{\tau }_0^2)\propto e^{-\frac{1}{2}{({\mu }_j-{\mu }_0)}^2/{\tau }_0^2}\hfill \\ \hfill p_0\left({\tau }_j^2\right|e_0,{\gamma }_0)& :\mathit{InvGa}(e_0,{\gamma }_0)\propto {\left({\tau }_j^2\right)}^{-(e_0+1)}·e^{-{\gamma }_0/{\tau }_j^2},p_0\left(\alpha \right|g_0,h_0):\mathit{Ga}(g_0,h_0)\propto {\alpha }^{(g_0-1)}·e^{-\alpha ·h_0}\hfill \end{array}$$

^*${\mathit{\beta }}_0,{\Sigma }_0\sim$ Gamma regression, ${\tilde{\mathit{\beta }}}_0,{\tilde{\Sigma }}_0\sim$ Logistic regression.

Appendix A.2. Posterior Inference for Outcome, Covariates, and Precision

Algorithm A1:Posterior inference ${\mathit{\theta }}_j^{*}=\{{\mathit{\beta }}_j^{*},{\sigma }_j^{2*},{\xi }_j^{*},{\tilde{\mathit{\beta }}}_j^{*}\}$

Require: initialize ${\mathit{\theta }}_j^{\left(old\right)}:\left\{\begin{array}{c}{\mathit{\beta }}_j\sim \mathit{MVN}({\mathit{\beta }}_0,{\sigma }_j^2{\Sigma }_0){\sigma }_j^2\sim \mathit{IG}(a_0,b_0){\xi }_j\sim \mathit{T}\left({\nu }_0\right){\tilde{\mathit{\beta }}}_j\sim \mathit{MVN}({\tilde{\mathit{\beta }}}_0,{\tilde{\Sigma }}_0)\hfill \end{array}\right.$ 1: repeat  2:     for $j=1,\cdots ,J$do▹ Assume J cluster memberships.   3:         Sample ${\mathit{\theta }}^{\left(new\right)}$ from the proposal densities $\mathit{q}$: ▹ Choose priors as $\mathit{q}$.               ${\mathit{\beta }}_j^{\left(new\right)}\sim {\mathit{q}}_{\beta }$, ${\sigma }_j^{2\left(new\right)}\sim {\mathit{q}}_{{\sigma }^2}$, ${\xi }_j^{\left(new\right)}\sim {\mathit{q}}_{\xi }$, ${\tilde{\mathit{\beta }}}_j^{\left(new\right)}\sim {\mathit{q}}_{\tilde{\beta }}$   4:         for ${\mathit{\theta }}_j^{\left(new\right)}$={${\mathit{\beta }}_j^{\left(new\right)},{\sigma }_j^{2\left(new\right)},{\xi }_j^{\left(new\right)},{\tilde{\mathit{\beta }}}_j^{\left(new\right)}$} do    5:            Compute the transition ratio, using the outcome models:                    $Rati{o}_{\theta }={\displaystyle \frac{{\prod }_{h=1}^{H}f\left(S_h\right|\mathit{X},{\mathit{\theta }}_j^{\left(new\right)})1·p_0\left({\mathit{\theta }}_j^{\left(new\right)}\right)·{\mathit{q}}_{\theta }\left({\mathit{\theta }}_j^{\left(old\right)}\right)}{{\prod }_{h=1}^{H}f\left(S_h\right|\mathit{X},{\mathit{\theta }}_j^{\left(old\right)})1·p_0\left({\mathit{\theta }}_j^{\left(old\right)}\right)·{\mathit{q}}_{\theta }\left({\mathit{\theta }}_j^{\left(new\right)}\right)}}$              Sample $U\sim \mathit{Unif}(0,1)$   6:            if $U<Rati{o}_{\theta }$then${\mathit{\theta }}_j^{*}={\mathit{\theta }}_j^{\left(new\right)}$ otherwise   ${\mathit{\theta }}_j^{*}={\mathit{\theta }}_j^{\left(old\right)}$   7:            end if  8:         end for  9:         Record ${\mathit{\theta }}_j^{*}$   10:     end for  11: until M posterior samples (${\mathit{\theta }}_{j=1,\cdots ,J}^{*}$) obtained. ▹ M is a sufficient sample size

¹ The outcome density is defined as: $f\left(S_h\right|\mathit{X},{\mathit{\theta }}_j)=\delta \left({\mathit{X}}_h^T{\tilde{\mathit{\beta }}}_j\right)\mathbb{1}(S_h=0)+\left[1-\delta \left({\mathit{X}}_h^T{\tilde{\mathit{\beta }}}_j\right)\right]f_{LSN}\left(S_h\right|{\mathit{X}}_h,{\mathit{\theta }}_j)$.

$$\begin{array}{cc}\hfill p\left({\pi }_j\right|c_0,& d_0,S,{\mathit{x}}_1):\mathit{Beta}(c_{new},d_{new})p\left({\mu }_j\right|{\mu }_0,{\tau }_0^2,S,{\mathit{x}}_2):\mathit{N}({\mu }_{new},{\tau }_{new}^2)\hfill \\ \hfill & \left\{\begin{array}{c}{c}_{new}=c_0+\sum _{h=1}^{n_j}{x}_{1h}\hfill \\ d_{new}=d_0+n_j-\sum _{h=1}^{n_j}{x}_{1h}\hfill \end{array}\right.\left\{\begin{array}{c}{\mu }_{new}=(n_j\overline{x_2}+{\mu }_0)/(n_j+1)\hfill \\ {\tau }_{new}^2={\tau }_j^2/(n_j+1)\hfill \end{array}\right.\hfill \\ \hfill p\left({\tau }_j^2\right|e_0,& {\gamma }_0,S,{\mathit{x}}_2):\mathit{InvGa}(e_{new},{\gamma }_{new})p\left(\alpha \right|g_0,h_0,h,J,\eta ,{\pi }_{\eta }):{\pi }_{\eta }\mathit{Ga}(g_0+J,h_0-log\left(\eta \right))\hfill \\ \hfill & +(1-{\pi }_{\eta })\mathit{Ga}(g_0+J-1,h_0-log\left(\eta \right))\hfill \\ & \left\{\begin{array}{c}{e}_{new}=e_0+n_j/2\hfill \\ {\gamma }_{new}={\gamma }_0+\frac{1}{2}\{\frac{n_j}{n_j+1}·{(\overline{x_2}-{\mu }_0)}^2+\sum _{h=1}^{n_j}{(x_{2h}-\overline{x_2})}^2\}\hfill \end{array}\right.\left\{\begin{array}{c}\eta |\alpha ,h\sim \mathit{Beta}(\alpha +1,h)\hfill \\ {\pi }_{\eta }={\displaystyle \frac{g_0+J-1}{g_0+J-1+h(h_0-log\left(\eta \right))}}\hfill \end{array}\right.\hfill \end{array}$$

Appendix C.1. Derivation of the Distribution of Precision α

In Section 4.1, the parameter model (posterior) of the precision term $\alpha$ is defined as

$$\begin{array}{cc}\hfill p\left(\alpha \right|J)\propto & p_0\left(\alpha \right)·{\alpha }^{J-1}·(\alpha +n)·\mathit{Beta}(\alpha +1,n)\hfill \\ \hfill p\left(\alpha \right|J,\eta ,g_0,h_0)\propto & {\pi }_{\eta }\mathit{Ga}(g_0+J,h_0-log\left(\eta \right))+(1-{\pi }_{\eta })\mathit{Ga}(g_0+J-1,h_0-log\left(\eta \right))\hfill \end{array}$$

To derive this, we first derive the distribution of the number of clusters given the precision parameter: $p\left(J\right|\alpha )$. Consider a trivial example where we want to determine the number of clusters that $n=5$ observations fall into. One possible arrangement would be that observations 1, 2, and 5 form new clusters, while observations 3 and 4 join an existing cluster. (note, the order is important):

• observation 1 forms a new cluster, probability = ${\displaystyle \frac{\alpha }{\alpha }}$
• observation 2 forms a new cluster, probability = ${\displaystyle \frac{\alpha }{\alpha +1}}$
• observation 3 enters into an existing cluster, probability = ${\displaystyle \frac{2}{\alpha +2}}$
• observation 4 enters into an existing cluster, probability = ${\displaystyle \frac{3}{\alpha +3}}$
• observation 5 forms a new cluster, probability = ${\displaystyle \frac{\alpha }{\alpha +4}}$

In this example, we have $J=3$ clusters. We want to find the probability of this arrangement. The probability is the following:

$$\begin{array}{cc}\hfill \left({\displaystyle \frac{\alpha }{\alpha }}\right)\left({\displaystyle \frac{\alpha }{\alpha +1}}\right)\left({\displaystyle \frac{2}{\alpha +2}}\right)\left({\displaystyle \frac{3}{\alpha +3}}\right)\left({\displaystyle \frac{\alpha }{\alpha +4}}\right)& \propto {\displaystyle \frac{{\alpha }^3}{\alpha (\alpha +1)(\alpha +2)(\alpha +3)(\alpha +4)}}\hfill \\ \hfill & ={\alpha }^3{\displaystyle \frac{\Gamma \left(\alpha \right)}{\Gamma (\alpha +5)}}\hfill \end{array}$$

Hence the probability of observing J clusters amongst a sample size of n is given by

$$p\left(J\right|\alpha )\propto {\alpha }^J{\displaystyle \frac{\Gamma \left(\alpha \right)}{\Gamma (\alpha +n)}}$$

This is also considered the likelihood function. The posterior on $\alpha$ is proportional to the likelihood times the prior, $p_0\left(\alpha \right)$.

$$\begin{array}{cc}\hfill p\left(\alpha \right|J)& \propto p\left(J\right|\alpha )p_0\left(\alpha \right)\hfill \\ \hfill & \propto {\alpha }^J{\displaystyle \frac{\Gamma \left(\alpha \right)}{\Gamma (\alpha +n)}}{p}_0\left(\alpha \right)\hfill \end{array}$$

The beta function $\mathit{Beta}(x,y)$ is defined as the following:

$$\begin{array}{cc}\hfill \mathit{Beta}(x,y)& ={\displaystyle \frac{\Gamma \left(x\right)\Gamma \left(y\right)}{\Gamma (x+y)}}\hfill \end{array}$$

We can find the beta function of $\alpha +1$ and n as follows:

$$\begin{array}{cc}\hfill \mathit{Beta}(\alpha +1,n)& ={\displaystyle \frac{\Gamma (\alpha +1)\Gamma \left(n\right)}{\Gamma (\alpha +1+n)}}\hfill \\ \hfill & \propto {\displaystyle \frac{\alpha \Gamma \left(\alpha \right)}{(\alpha +n)\Gamma (\alpha +n)}}\hfill \\ \hfill {\displaystyle \frac{\Gamma \left(\alpha \right)}{\Gamma (\alpha +n)}}& \propto \mathit{Beta}(\alpha +1,n){\displaystyle \frac{\alpha +n}{\alpha }}\hfill \end{array}$$

Thus the posterior simplifies to the following:

$$\begin{array}{cc}\hfill p\left(\alpha \right|J)& \propto {\alpha }^J·\mathit{Beta}(\alpha +1,n)·{\displaystyle \frac{\alpha +n}{\alpha }}·p_0\left(\alpha \right)\hfill \\ \hfill & \propto p_0\left(\alpha \right)·{\alpha }^{J-1}·(\alpha +n)·\mathit{Beta}(\alpha +1,n)\hfill \end{array}$$

Now, under the $\mathit{Ga}(g_0,h_0)$ prior for $\alpha$, substituting $p_0\left(\alpha \right)$ with $\mathit{Ga}(g_0,h_0)$, then

$$\begin{array}{cc}\hfill p\left(\alpha \right|J,\eta ,g_0,h_0)\propto & {\alpha }^{g_0+j-2}·(\alpha +n)·e^{-\alpha (h_0-log\left(\eta \right)}\hfill \\ \hfill \propto & {\pi }_{\eta }\mathit{Ga}(g_0+J,h_0-log\left(\eta \right))+(1-{\pi }_{\eta })\mathit{Ga}(g_0+J-1,h_0-log\left(\eta \right))\hfill \end{array}$$

Appendix C.2. Outcome Data Model of S h Development with MAR Covariate x 1 for the Discrete Clusters

Prior to the outcome parameter estimation, the missing covariates should be imputed first to obtain the complete covariate model beforehand. In this study, if the binary covariate $x_{1h}$ is the only covariate with missingness, we develop the imputation model to impute the binary covariate $x_{1h}$, taking the following steps below, then update $\mathit{\beta },{\sigma }^2,\xi ,\tilde{\mathit{\beta }}$ based on the posterior sampling detailed in Algorithm (A1) in Appendix (A.2). The imputation model for $x_{1h}$ is approximated by the joint:

$$f\left(x_{1h}\right|S_h,x_{2h},{\mathit{\beta }}_j,{\sigma }_j,{\xi }_j,\widetilde{{\mathit{\beta }}_j},{\pi }_j)\propto f(S_h,x_{1h}|x_{2h},{\mathit{\beta }}_j,{\sigma }_j,{\xi }_j,\widetilde{{\mathit{\beta }}_j},{\pi }_j)$$

where

$$\begin{array}{cc}\hfill f(S_h,x_{1h}|& x_{2h},{\mathit{\beta }}_j,{\sigma }_j^2,{\xi }_j,\widetilde{{\mathit{\beta }}_j},{\pi }_j)=f\left(S_h\right|x_{1h},x_{2h},{\mathit{\beta }}_j,{\sigma }_j^2,{\xi }_j,\widetilde{{\mathit{\beta }}_j})·f_{Bern}\left(x_{1h}\right|{\pi }_j)\hfill \\ \hfill & =\delta \left({\mathit{X}}_h^T{\tilde{\mathit{\beta }}}_j\right)\mathbb{1}(S_h=0)·{\pi }_j^{x_{1h}}{\left(1-{\pi }_j\right)}^{1-x_{1h}}+\left[1-\delta \left({\mathit{X}}_h^T{\tilde{\mathit{\beta }}}_j\right)\right]{\displaystyle \frac{2}{{\sigma }_j{S}_h}}\hfill \\ \hfill & ·\varphi \left(\frac{log{S}_h-{\mathit{X}}^T{\mathit{\beta }}_j}{{\sigma }_j}\right)·\mathrm{\Phi }\left({\xi }_j\frac{log{S}_h-{\mathit{X}}^T{\mathit{\beta }}_j}{{\sigma }_j}\right)·{\pi }_j^{x_{1h}}{\left(1-{\pi }_j\right)}^{1-x_{1h}}\hfill \end{array}$$

which serves as the joint density that we can use to sample the imputation values. For example,

$$\begin{array}{cc}\hfill f_{Bern}& (x_{1h}=1|S_h,x_{2h},{\mathit{\beta }}_j,{\sigma }_j^2,{\xi }_j,\widetilde{{\mathit{\beta }}_j},{\pi }_j)\propto f(S_h,x_{1h}=1|x_{2h},{\mathit{\beta }}_j,{\sigma }_j^2,{\xi }_j,\widetilde{{\mathit{\beta }}_j},{\pi }_j)\hfill \\ \hfill & =\delta ({\tilde{\beta }}_{j0}+{\tilde{\beta }}_{j1}+{\tilde{\beta }}_{j2}{x}_{2h})\mathbb{1}(S_h=0)·{\pi }_j+\left[1-\delta ({\tilde{\beta }}_{j0}+{\tilde{\beta }}_{j1}+{\tilde{\beta }}_{j2}{x}_{2h})\right]{\displaystyle \frac{2}{{\sigma }_j{S}_h}}\hfill \\ \hfill & ·\varphi \left(\frac{log{S}_h-({\beta }_{j0}+{\beta }_{j1}+{\beta }_{j2}{x}_{2h})}{{\sigma }_j}\right)·\mathrm{\Phi }\left({\xi }_j\frac{log{S}_h-({\beta }_{j0}+{\beta }_{j1}+{\beta }_{j2}{x}_{2h})}{{\sigma }_j}\right){\pi }_j\hfill \\ \hfill f_{Bern}& (x_{1h}=0|S_h,x_{2h},{\mathit{\beta }}_j,{\sigma }_j^2,{\xi }_j,\widetilde{{\mathit{\beta }}_j},{\pi }_j)\propto f(S_h,x_{1h}=0|x_{2h},{\mathit{\beta }}_j,{\sigma }_j^2,{\xi }_j,\widetilde{{\mathit{\beta }}_j},{\pi }_j)\hfill \\ \hfill & =\delta ({\tilde{\beta }}_{j0}+{\tilde{\beta }}_{j2}{x}_{2h})\mathbb{1}(S_h=0)·(1-{\pi }_j)+\left[1-\delta ({\tilde{\beta }}_{j0}+{\tilde{\beta }}_{j2}{x}_{2h})\right]{\displaystyle \frac{2}{{\sigma }_j{S}_h}}\hfill \\ \hfill & ·\varphi \left(\frac{log{S}_h-({\beta }_{j0}+{\beta }_{j2}{x}_{2h})}{{\sigma }_j}\right)·\mathrm{\Phi }\left({\xi }_j\frac{log{S}_h-({\beta }_{j0}+{\beta }_{j2}{x}_{2h})}{{\sigma }_j}\right)·(1-{\pi }_j)\hfill \end{array}$$

Then, we can impute $x_{1h}$ with the values sampled from $Bern\left({\pi }_{x_1}^{*}\right)$ where

$${\pi }_{x_1}^{*}=\frac{f(S_h,x_{1h}=1|x_{2h},{\mathit{\beta }}_j,{\sigma }_j^2,{\xi }_j,\widetilde{{\mathit{\beta }}_j},{\pi }_j)}{f(S_h,x_{1h}=1|x_{2h},{\mathit{\beta }}_j,{\sigma }_j^2,{\xi }_j,\widetilde{{\mathit{\beta }}_j},{\pi }_j)+f(S_h,x_{1h}=0|x_{2h},{\mathit{\beta }}_j,{\sigma }_j^2,{\xi }_j,\widetilde{{\mathit{\beta }}_j},{\pi }_j)}$$

Note that in R, the computation can be difficult when the numerator is too small. We suggest the following tricks.

$$\begin{array}{cc}\hfill p_1=f(S_h,x_{1h}=1|& x_{2h},{\mathit{\beta }}_j,{\sigma }_j^2,{\xi }_j,\widetilde{{\mathit{\beta }}_j},{\pi }_j)\hfill \\ \hfill p_0=f(S_h,x_{1h}=0|& x_{2h},{\mathit{\beta }}_j,{\sigma }_j^2,{\xi }_j,\widetilde{{\mathit{\beta }}_j},{\pi }_j)\hfill \\ \hfill & {\pi }_{x_1}^{*}={\displaystyle \frac{e^{log\left(p_1\right)}}{e^{log\left(p_1\right)}+e^{log\left(p_0\right)}}}·{\displaystyle \frac{e^{-log\left(p_1\right)}}{e^{-log\left(p_1\right)}}}={\displaystyle \frac{1}{1+e^{log\left(p_0\right)-log\left(p_1\right)}}}\hfill \end{array}$$

Finally, the outcome model that is required to compute the parameter $\mathit{\theta }=\{{\mathit{\beta }}_j,{\sigma }_j^2,{\xi }_j,\widetilde{{\mathit{\beta }}_j}\}$ in Metropolis-Hastings in Algorithm (A1) is obtained by summing the joint of $S_h$ and $x_{1h}$ (marginalize) out the MAR covariate $x_{1h}$, shown in Equation (9), as below.

$$\begin{array}{cc}\hfill f\left(S_h\right|x_{2h},& {\mathit{\beta }}_j,{\sigma }_j^2,{\xi }_j,\widetilde{{\mathit{\beta }}_j},{\pi }_j)=\sum _{x_{1h}=0}^{1}f(S_h,x_{1h}|x_{2h},{\mathit{\beta }}_j,{\sigma }_j^2,{\xi }_j,\widetilde{{\mathit{\beta }}_j},{\pi }_j)\hfill \\ \hfill & =f(S_h,x_{1h}=1|x_{2h},{\mathit{\beta }}_j,{\sigma }_j^2,{\xi }_j,\widetilde{{\mathit{\beta }}_j},{\pi }_j)+f(S_h,x_{1h}=0|x_{2h},{\mathit{\beta }}_j,{\sigma }_j^2,{\xi }_j,\widetilde{{\mathit{\beta }}_j},{\pi }_j)\hfill \end{array}$$

Appendix C.3. Covariate Data Model of x 2 Development with MAR Covariate x 1 for the Continuous Clusters

The parameter-free distributions $f_0\left(y\right|x)$ and $f_0\left(x\right)$ as data models for continuous clusters are needed to calculate the probabilities of cluster membership and for the post-processing calculations for prediction in the DPM. However, when MAR covariates are present, it gives extra complexity in specifying distribution to integrate out the parameters. Recall the integrals we are attempting to find are the following:

$$\begin{array}{c}\hfill f_0\left(x_i\right)=\int f\left(x_i\right|\mathit{w})d{G}_0\left(\mathit{w}\right)=\int f\left(x_i\right|\mathit{w})p\left(\mathit{w}\right)d\mathit{w}\end{array}$$

If binary covariate $x_1$ is missing, then we will need to replace the distribution $f\left(x\right|\mathit{w})$ with the continuous distribution (Gaussian) of $x_2$, which is $f\left(x_2\right|{\mu }_j,{\tau }_j^2)$. The derivation of the parameter-free distribution $f_0\left(x_1\right)$ and $f_0\left(x_2\right)$ for the continuous cluster is as below.

$$\begin{array}{cc}\hfill f_0\left(x_1\right)& \\ \hfill & =\int f\left(x_1\right|\pi )p\left(\pi \right)d\mu d\pi \hfill \\ \hfill & =\int {\pi }^{x_1}{\left(1-\pi \right)}^{1-x_1}{\displaystyle \frac{1}{\mathit{Beta}(c_0,d_0)}}{\pi }^{(c_0-1)}{(1-\pi )}^{(d_0-1)}d\pi \hfill \\ \hfill & ={\displaystyle \frac{1}{\mathit{Beta}(c_0,d_0)}}\int {\pi }^{(x_1+c_0-1)}{(1-\pi )}^{(1-x_1+d_0-1)}d\pi \hfill \\ \hfill & =\frac{\mathit{Beta}(x_1+c_0,1-x_1+d_0)}{\mathit{Beta}(c_0,d_0)}·\underset{=1,\mathrm{beta}\mathrm{distribution}}{\underbrace{\int {\displaystyle \frac{{\pi }^{(x_1+c_0-1)}{(1-\pi )}^{(1-x_1+d_0-1)}}{\mathit{Beta}(x_1+c_0,1-x_1+d_0)}}d\pi }}\hfill \end{array}$$

$$\begin{array}{cc}\hfill f_0& \left(x_2\right)\hfill \\ \hfill & =\int \int f\left(x_2\right|\mu ,{\tau }^2)p\left(\mu \right|{\tau }^2)p\left({\tau }^2\right)d\mu d{\tau }^2\hfill \\ \hfill & =\int \int {\displaystyle \frac{1}{\sqrt{2\pi {\tau }^2}}}exp\left\{-{\displaystyle \frac{1}{2{\tau }^2}}{\left(x_2-\mu \right)}^2\right\}\times {\displaystyle \frac{1}{\sqrt{2\pi {\tau }^2}}}exp\left\{-{\displaystyle \frac{1}{2{\tau }^2}}{\left(\mu -{\mu }_0\right)}^2\right\}\hfill \\ \hfill & \times {\displaystyle \frac{{\gamma }_0^{e_0}}{\Gamma \left(e_0\right)}}{\left({\tau }^2\right)}^{-e_0-1}{e}^{-{\gamma }_0/{\tau }^2}d\mu d{\tau }^2\hfill \\ \hfill & ={\displaystyle \frac{{\gamma }_0^{e_0}}{2\pi \Gamma \left(e_0\right)}}\int \int {\left({\tau }^2\right)}^{-e_0-2}exp\left\{-{\displaystyle \frac{1}{2{\tau }^2}}{\left(x_2-\mu \right)}^2-{\displaystyle \frac{1}{2{\tau }^2}}{\left(\mu -{\mu }_0\right)}^2-{\displaystyle \frac{{\gamma }_0}{{\tau }^2}}\right\}d\mu d{\tau }^2\hfill \end{array}$$

The first step is to integrate with respect to $\mu$. First, we’ll simplify the exponent.

$$\begin{array}{cc}\hfill -{\displaystyle \frac{1}{2{\tau }^2}}{\left(x_2-\mu \right)}^2& -{\displaystyle \frac{1}{2{\tau }^2}}{\left(\mu -{\mu }_0\right)}^2-{\displaystyle \frac{{\gamma }_0}{{\tau }^2}}\hfill \\ \hfill & =-{\displaystyle \frac{1}{2{\tau }^2}}\left[x_2^2-2x_2\mu +{\mu }^2+{\mu }^2-2{\mu }_0\mu +{\mu }_0^2\right]-{\displaystyle \frac{{\gamma }_0}{{\tau }^2}}\hfill \\ \hfill & =-{\displaystyle \frac{1}{2{\tau }^2}}\left[2{\mu }^2-2(x_2+{\mu }_0)\mu \right]-{\displaystyle \frac{1}{2{\tau }^2}}\left[x_2^2+{\mu }_0^2\right]-{\displaystyle \frac{{\gamma }_0}{{\tau }^2}}\hfill \\ \hfill & =-{\displaystyle \frac{2}{2{\tau }^2}}\left[{\mu }^2-(x_2+{\mu }_0)\mu +{\displaystyle \frac{{(x_2+{\mu }_0)}^2}{4}}\right]+{\displaystyle \frac{1}{{\tau }^2}}\left({\displaystyle \frac{{(x_2+{\mu }_0)}^2}{4}}\right)\hfill \\ \hfill & -{\displaystyle \frac{x_2^2+{\mu }_0^2}{2{\tau }^2}}-{\displaystyle \frac{{\gamma }_0}{{\tau }^2}}\hfill \\ \hfill & =-{\displaystyle \frac{1}{2({\tau }^2/2)}}{\left(\mu -{\displaystyle \frac{x_2+{\mu }_0}{2}}\right)}^2+{\displaystyle \frac{{(x_2+{\mu }_0)}^2}{4{\tau }^2}}-{\displaystyle \frac{x_2^2+{\mu }_0^2}{2{\tau }^2}}-{\displaystyle \frac{{\gamma }_0}{{\tau }^2}}\hfill \end{array}$$

The integrand will have the kernel of a normal distribution for $\mu$ with mean ${\displaystyle \frac{x_2+{\mu }_0}{2}}$ and variance ${\displaystyle \frac{{\tau }^2}{2}}$.

$$\begin{array}{cc}\hfill f_0& \left(x_2\right)\hfill \\ \hfill & ={\displaystyle \frac{{\gamma }_0^{e_0}}{2\pi \Gamma \left(e_0\right)}}\int \underset{\mathrm{term}\mathrm{from}\mu \mathrm{integral}}{\underbrace{\sqrt{2\pi ({\tau }^2/2)}}}\times {\left({\tau }^2\right)}^{-e_0-2}\times exp\left\{{\displaystyle \frac{{(x_2+{\mu }_0)}^2}{4{\tau }^2}}-{\displaystyle \frac{x_2^2+{\mu }_0^2}{2{\tau }^2}}-{\displaystyle \frac{{\gamma }_0}{{\tau }^2}}\right\}d{\tau }^2\hfill \\ \hfill & ={\displaystyle \frac{{\gamma }_0^{e_0}}{2\sqrt{\pi }\Gamma \left(e_0\right)}}\int {\left({\tau }^2\right)}^{-e_0-3/2}exp\left\{-{\displaystyle \frac{1}{{\tau }^2}}\left(-{\displaystyle \frac{x_2^2+2x_2{\mu }_0+{\mu }_0^2}{4}}+{\displaystyle \frac{x_2^2+{\mu }_0^2}{2}}+{\gamma }_0\right)\right\}d{\tau }^2\hfill \\ \hfill & ={\displaystyle \frac{{\gamma }_0^{e_0}}{2\sqrt{\pi }\Gamma \left(e_0\right)}}\int {\left({\tau }^2\right)}^{-e_0-1/2-1}exp\left\{-{\displaystyle \frac{1}{{\tau }^2}}\left({\displaystyle \frac{{(x_2^2-{\mu }_0)}^2}{4}}+{\gamma }_0\right)\right\}d{\tau }^2\hfill \end{array}$$

The integrand is the kernel of an inverse gamma distribution with shape parameter $e_0+{\displaystyle \frac{1}{2}}$ and scale parameter ${\displaystyle \frac{{(x_2^2-{\mu }_0)}^2}{4}}+{\gamma }_0$.

$$\begin{array}{cc}\hfill f_0\left(x_2\right)& ={\displaystyle \frac{{\gamma }_0^{e_0}}{2\sqrt{\pi }\Gamma \left(e_0\right)}}\times \Gamma (e_0+1/2){\left({\displaystyle \frac{{(x_2^2-{\mu }_0)}^2}{4}}+{\gamma }_0\right)}^{-e_0-1/2}\hfill \end{array}$$

As shown above, a closed-form expression can be determined, but it is not always the case since it can be extremely complicated. To simplify, we instead might have to consider a Monte Carlo integral.