# Analyzing covariate clustering effects in healthcare cost subgroups:
insights and applications for prediction

Zhengxiao Li 111School of Insurance and Economics, University of International
Business and Economics, Beijing, China. Yifan Huang 222corresponding author:
School of Insurance and Economics, University of International Business and
Economics, Beijing, China. Email<EMAIL_ADDRESS>Yang Cao 333School of
Insurance and Economics, University of International Business and Economics,
Beijing, China.

###### Abstract

Healthcare cost prediction is a challenging task due to the high-
dimensionality and high correlation among covariates. Additionally, the
skewed, heavy-tailed, and often multi-modal nature of cost data can complicate
matters further due to unobserved heterogeneity. In this study, we propose a
novel framework for finite mixture regression models that incorporates
covariate clustering methods to better account for the effects of clustered
covariates on subgroups of the outcome, which enables a more accurate
characterization of the complex distribution of the data. The proposed
framework can be formulated as a convex optimization problem with an
additional penalty term based on the prior similarity of the covariates. To
efficiently solve this optimization problem, a specialized EM-ADMM algorithm
is proposed that integrates the alternating direction multiplicative method
(ADMM) into the iterative process of the expectation-maximizing (EM)
algorithm. The convergence of the algorithm and the efficiency of the
covariate clustering method are verified using simulation data, and the
superiority of the approach over traditional regression techniques is
demonstrated using two real Chinese medical expenditure datasets. Our
empirical results provide valuable insights into the complex network graph of
the covariates and can inform business practices, such as the design and
pricing of medical insurance products.
Keywords: covariate clustering; finite mixture regression models; EM-ADMM
algorithm; healthcare costs prediction.

Article History: August 28, 2024

##  1 Introduction

The prediction of healthcare costs for individuals is a crucial risk
management concern for healthcare systems and health insurers in determining
appropriate premiums (Fellingham et al., 2015, Xie et al., 2016). While
regression methods have been well-developed for predictive modelling of
healthcare costs from individual clinical data in the economic and biomedical
statistics literature (Bertsimas et al., 2008, Neelon et al., 2015, Cao et
al., 2012, Duncan et al., 2016), high-dimensional regression models are often
necessary due to the large number of risk factors in healthcare datasets, some
of which are highly correlated and grouped into sets (or clusters) known to
the data analyst (e.g. demographic information such as age, gender, income;
health status such as smoking, drinking, surgery history; and clinical
measurements such as blood pressure, blood glucose, resting heart rate).
Despite successful covariate selection, the covariates need further clustering
to increase model accuracy and interpretability (Andrade et al., 2021).
Another challenge is that covariates respond differently to the response
variable due to unobserved heterogeneity and latent factors, and common
distributions such as Gamma, Weibull, and Lognormal may not fully capture the
heavy-tailedness and multi-modality of individual healthcare costs data
(Mihaylova et al., 2011), which also suggests the need to identify the
subgroups of the outcome data. To address these challenges, this paper
proposes a statistical regression model for healthcare cost prediction, with
the aim of precise covariate selection and covariate clustering for
heterogeneous subgroups of individuals.

Given the heterogeneity which characterizes individual clinical data, many
approaches have been proposed in the statistical literature to identify
unobserved subgroups from a heterogeneous population (Ma and Huang, 2017). One
such approach is to treat the outcome data as a mixture of several subgroups
with their own parameters and to fit a finite mixture regression model (FMR).
Previous studies have explored FMR models with no covariates (Ickowicz and
Sparks, 2017, Kurz and Hatfield, 2019), as well as those with low-dimensional
covariates (Shen and He, 2015). While mixture models can easily incorporate
covariate effects due to their theoretical framework, they are limited when
covariates are highly correlated in high-dimensional problems. To address this
issue, regularization approaches have been developed in FMR models (Khalili
and Chen, 2007, Khalili and Lin, 2013), with the $\ell_{p}$ regularization
term being a particularly popular choice. More recently, a high-dimensional
mixture Gamma regression model was proposed to fit heavy-tailed insurance
claim data with neural networks, allowing for covariate selection and
detection of non-linear regression patterns (Devriendt et al., 2021). This
indicates a mixture Gamma distribution is preferred to model the right-skewed,
heavy-tailed and multi-modal insurance claim data (Yin and Lin, 2016,
Miljkovic and Grün, 2016). However, these methods cannot reveal the structure
of the relationship between highly correlated covariates, although they can
reduce the dimensionality of the data. Therefore, further clustering methods
for highly correlated covariates are still needed to improve the
interpretability of the model.

One commonly used approach to reducing the dimensionality of covariates is
unsupervised clustering. For instance, Tsai and Cheng (2021) applied four
clustering approaches (i.e., Ward’s hierarchical clustering, divisive
hierarchical clustering, K-means clustering, and Gaussian mixture model
clustering) to improve forecasting performance of mortality models in
insurance aspect. The drawback of these methods is that it is generally
performed independently of regression, overlooking the response-covariate
relationships, leading to undesirable clustering results. A approach for
identifying clusters among the covariates and simultaneously estimating the
regression coefficients has been proposed in several studies, such as Witten
et al. (2014), Sui et al. (2018), Hallac et al. (2015) and Andrade et al.
(2021). The key idea of Andrade et al. (2021) is that while taking into
account the response information is to formulate a convex classification
problem with a pair-wise group lasso penalty, which selectively shrinks the
coefficients of such covariates that are highly correlated but not related to
the origin. This modelling framework can be also referred to a network lasso
(Hallac et al., 2015) or generalized lasso problem (Zhu, 2017). However, the
existing methods have limitations in identifying and interpreting different
clustering results from various subgroups of healthcare costs data.

In this paper, we present a new statistical framework for healthcare cost
modeling that identifies the contribution of clustered covariates to the
subgroups of the outcome, with clustering results varying across different
components of the FMR model. The proposed approach builds on the work of
Andrade et al. (2021), which clusters the regression coefficients of
covariates using pairwise fused penalty. Specifically, the framework starts
with an FMR model based on the exponential distribution family, where the
outcome is treated as a mixture of different subgroups. Each subgroup is
assumed to be an underlying distribution with its own set of parameters. The
mean parameter of each underlying distribution is formulated by linearly
combining the covariates using a log-link function. To obtain the covariate
clustering results during mixture regression, prior similarity information
about the covariates is used, and two convex penalty terms are added to the
log-likelihood function. These penalty terms account for the fit deviation of
observed samples (the $\ell_{2}$ norm of the coefficients) and the distance of
similarity between covariates (the $\ell_{1}$ norm of the difference of pair-
wise coefficients) in different subgroups. The proposed model adaptively
explores the pair-wise structure of a large number of highly correlated
covariates by shrinking the group-specific coefficients in different mixture
components. The approach can cluster the covariates information when subgroup
regression is performed simultaneously, which is an improvement over
traditional methods that separate the two steps and model them disjointly. The
proposed framework can also be applied to one-dimensional or non-mixture
models when there is weak heterogeneity within the observed sample. Notably,
the proposed model can be viewed as an extension of the linear regression
models with convex clustering method proposed by Cheng et al. (2022).

Our second contribution is employing a specialized EM-ADMM algorithm for
optimizing the complicated objective function. In this algorithm, we perform
an alternating direction method of multipliers (ADMM) at the M-step of the
Expectation-Maximization (EM) algorithm, instead of using the Newton-Raphson
optimization procedure to estimate the unknown parameters. The ADMM algorithm
has good convergence properties for convex objective functions with the
$\ell_{p}\ (p\geq 1)$ regularization terms, as demonstrated by Boyd et al.
(2011) and Chi and Lange (2015). Furthermore, the $\ell_{1}$ regularization
can shrink some pairwise differences of the parameter estimates to zero. To
obtain satisfactory covariate clustering results for each component in the
mixture model, we choose the penalty imposed on the regression coefficients
within each component of the FMR model to be proportional to the mixture
weight, following the work of Khalili and Chen (2007). Our simulation study
shows that the proposed algorithm provides a parallel computational strategy
for the parameter estimation of each component in the mixture model, resulting
in a fast convergence rate and high accuracy.

We present two practical healthcare costs datasets to demonstrate the
effectiveness of our proposed model. The first dataset pertains to a health
insurance product developed for chronic diabetes, with the primary goal being
the prediction of individual medical costs associated with diabetes based on
personal patient information. This dataset is unique in that it contains a
questionnaire consisting of 111 questions, many of which are similar and can
be grouped into specific sets of highly correlated covariates measuring self-
perception of diabetes for each patient. Given that medical utilization costs
for diabetes patients are strongly correlated with comorbidities and
complications, patients with a better understanding of diabetes are more
likely to receive timely treatment, including regular glucose-lowering
medication or insulin injections, thus reducing the probability of
complications and related costs. Our motivation to apply the covariate
clustering method is to uncover the relationship structures within the
questionnaire data, which can help improve the model’s predictive performance.
The second dataset is short-term medical insurance claims data provided by a
Chinese insurance company, where the ICD codes can be clustered into specific
groups for different subgroups during mixture regression. Our empirical
results demonstrate that our proposed model is suitable for heavy-tailed and
multi-modal medical costs data, and that the combination of regression and
covariate clustering can significantly improve model accuracy and
interpretability for predicting future insurance claims. This can be
particularly useful in guiding business practices such as designing and
pricing medical insurance products.

The rest of this paper is structured as follows: in Section 2, we first define
the framework of FMR models with covariate clustering, including the symbols
and main objective function. An EM-ADMM algorithm for solving the optimization
problem as well as some computation details are then shown in Section 3.
Section 4 verifies the performance of proposed model through simulation
studies, while in Section 5 empirical analyses of two real datasets are
illustrated. Conclusions are given in Section 6. To make the framework and
algorithm more accessible to practitioners, we provide the R programming codes
for implementation at https://github.com/huangyf2217/FMR-covariates-
clustering.

##  2 Finite mixture of regression model with covariate clustering

###  2.1 Finite mixture of regression models

Finite mixture models provide a tool for analyzing data suspected of coming
from a heterogeneous population. The positive continuous response variable of
interest (e.g., healthcare costs) is consistently denoted as $Y$, and
$\bm{x}=(x_{1},x_{2},...,x_{p})$ is the $p$-dimensional covariate vector
related to the risk.

Generally, the unobservable latent factors in heterogeneous response variable
can be distinguished by the mixture of underlying distributions. We first
assume that the random variable $Y$ comes from a parametric exponential family
$\mathbb{F}=\left\\{f(y;\mu,\phi);(\mu,\phi)\in\mu\times\left(0,+\infty\right)\right\\}$,
where $\mu$ is the specific mean parameter and $\phi$ is the dispersion
parameter. In generalized linear models (GLMs), the mean parameter can be
expressed as a function related to the covariates, specifically,
$\mu:=\mu(\bm{x})$. Then, the conditional density function of $Y$ given
$\bm{x}$ has the form, extending the above structure to an FMR model:

$f_{\text{mix}}\left(y\left|\bm{x},\bm{\omega},\bm{\beta}_{0},\mathbf{B},\bm{\phi}\right.\right)=\sum_{h=1}^{H}{\omega_{h}}{f\left(y\left|\mu_{h}(\bm{x}),\phi_{h}\right.\right)},$
(2.1)

where $\mu_{h}(\bm{x})=g^{-1}(\beta_{0,h}+\bm{x}^{T}\bm{\beta}_{h})$ with a
link function $g(\cdot)$.
$\bm{\beta}_{0}=\left(\beta_{0,1},\beta_{0,2},\dots,\beta_{0,H}\right)^{T}\in\mathbb{R}^{1\times
H}$ and
$\mathbf{B}=(\bm{\beta}_{1},\bm{\beta}_{2},\dots,\bm{\beta}_{H})^{T}\in\mathbb{R}^{p\times
H}$ contain the intercept terms and regression coefficients respectively,
where $\bm{\beta}_{h}=(\beta_{1,h},\beta_{2,h},\dots,\beta_{p,h})$ is the
vector of regression coefficients in the $h$ component.
$\bm{\phi}=\left(\phi_{1},\phi_{2},\dots,\phi_{H}\right)^{T}\in\mathbb{R}^{1\times
H}$ represents the dispersion parameters of all $H$ components. The component
weight $\bm{\omega}=\left(\omega_{1},\omega_{2},\dots,\omega_{H}\right)^{T}$
satisfies the constraints $\omega_{h}>0$ and $\sum_{h=1}^{H}\omega_{h}=1$, so
the mixture density function is actually the weighted sum of all components’
densities. The form of $f(\cdot)$ can vary according to the response variable.
Specifically, a Gamma distribution is often chosen as the underlying
distribution in FMR for fitting the positive continuous data in insurance
applications.

Given a sample of observations
$(\bm{x}_{1},y_{1}),(\bm{x}_{2},y_{2}),...,(\bm{x}_{n},y_{n})$, the
(penalized) log-likelihood function for the FMR model is as follows:

$\breve{\ell}_{n}(\bm{\Psi})=\ell_{n}(\bm{\Psi})-p_{n}(\bm{\Psi}),$ (2.2)

and

$\ell_{n}(\bm{\Psi})=\sum_{i=1}^{n}\log
f_{\text{mix}}\left(y_{i}\left|\bm{x}_{i},\bm{\omega},\bm{\beta}_{0},\mathbf{B},\bm{\phi}\right.\right)=\sum_{i=1}^{n}\log\sum_{h=1}^{H}{\omega_{h}}{f\left(y_{i}\left|\mu_{h}(\bm{x}_{i}),\phi_{h}\right.\right)},$
(2.3)

where
$\bm{\Psi}:=\left(\bm{\omega},\bm{\beta}_{0},\bm{\beta},\bm{\phi}\right)$,
$x_{i}=(x_{1,p},...,x_{i,p})$ is the $p$-dimensional covariates vector, and
the regularization term is defined as
$p_{n}(\bm{\Psi})=\gamma\sum_{h=1}^{H}{\omega_{h}}{\left\|\bm{\beta}_{h}\right\|_{2}^{2}}$
with $\gamma$ being a non-negative regularization parameter. This is in line
with the common practice of linking the penalty term to sample size, and
improves the model’s ability to identify important covariates, given that the
weight $\omega_{h}$ reflects to some extent the size of the $h$th
subpopulation (Khalili and Chen, 2007). The weighted $\ell_{2}$ norm is able
to compress the insignificant parameters in each component to avoid the model
overfitting. In addition to quantifying the probability that the observation
belongs to each component, the regularized FMR models are able to describe how
the covariates are related to the response variable under the influence of
various latent factors. The number of subgroups $H$ and the corresponding
underlying distributions are usually pre-specified, which is usually not too
large in healthcare costs modelling.

###  2.2 The framework of convex covariate clustering

For analysis of the covariate clustering effects on different components in
the FMR model, we can identify the same (or similar) values of corresponding
regression coefficients in each component, thus grouping the “similar”
covariates with similar semantics or values.

Suppose a similarity matrix $\mathbf{S}=\\{s_{j,k}\\}_{p\times p}$ is given a
priori, where the element $s_{j,k}$ defines the similarity measure between the
covariate $j$ and $k$ that are known to the experts. Generally,
$s_{j,k}=s_{k,j}$ is held. In the graphical interpretation of clustering, the
similarity matrix can be regarded as an undirected weighted graph, in which
each covariate corresponds to a node in a graph, and an edge connects nodes
$j$ and $k$ whenever $s_{j,k}>0$.

Following the idea of Andrade et al. (2021), we introduce an additional
similarity regularization term in (2.2), where the penalized objective
function is defined as

$\widetilde{\ell}_{n}(\bm{\Psi})=\breve{\ell}_{n}(\bm{\Psi})-\frac{v}{2}\sum_{h=1}^{H}\sum_{j,k\in
I}\omega_{h}s_{j,k}\lvert{\beta}_{j,h}-{\beta}_{k,h}\rvert,$ (2.4)

where $v$ is a non-negative tuning constant, $I=\left\\{1,2,...,p\right\\}$ is
the index set of covariates, $\lvert{\beta}_{j,h}-{\beta}_{k,h}\rvert$
measures the Euclidean distance between any pair of two covariates $j$ and $k$
that both have effects on the $h$th component of the outcome. The similarity
penalty is weighted by the value of $s_{j,k}$. We can cluster two covariates
$j$ and $k$, such that $\lvert\hat{\beta}_{j,h}-\hat{\beta}_{k,h}\rvert$ is
arbitrarily small when $v$ is sufficiently large, where
$\hat{\bm{\beta}_{h}}=(\hat{\beta}_{1,h},...,\hat{\beta}_{p,h})$ is the
maximization of (2.4). Note that the two regularization terms imposed on the
regression coefficients both are set to be proportional to the component
weight $\omega_{h}$ in (2.4), resulting in the component-specific clustering
result of the covariates. By simply setting $H=1$ in the objective function
above, the proposed framework can also be easily regarded as the traditional
GLMs with covariate clustering method.

##  3 Parameter estimation method

###  3.1 The optimization process of EM-ADMM algorithm

The two regularization terms in the penalized objective function ensure the
identifiability of our model. When $\gamma$ and $v$ are greater than zero, the
optimization problem (2.4) is strongly convex, thus we are guaranteed to find
the global optima by using appropriate algorithms.

For simple FMR models, the Expectation-Maximization (EM) algorithm provides a
convenient tool to solve the problem of the subgroup analysis based on the
maximum likelihood estimation method (Dempster et al., 1977). This algorithm
is generally divided into two steps: First, computing the conditional
expectation of the complete log-likelihood function for a given sample data
belonging to the specific latent subgroup (i.e. the E-step); second,
maximizing the conditional expectation of complete log-likelihood and
iteratively updating the estimates of $\bm{\Psi}$ (i.e. the M-step). However,
the commonly used Newton-type method is inadequate for solving the
optimization problem with multiple pairwise penalty terms, whereas the
difficulty of M-step optimization depends on the form of the expected complete
log-likelihood.

In this section, we discuss a numerical optimization method that uses the EM
algorithm to estimate the parameters in the FMR model, but incorporates the
Alternating Direction Method of Multipliers (ADMM) algorithm in the M-step to
solve the maximization problem. We call it EM-ADMM algorithm. Consistent with
the EM algorithm, we first introduce the latent variables $u_{i,h}$ into the
objective function (2.4), where $u_{i,h}$ indicates the class of the $i$th
observation, and a certain $u_{i,h}=1$ if and only if $y_{i}$ comes from the
$h$th mixture component. Then the penalized complete log-likelihood function
based on the latent variables takes the following form:

$\displaystyle\widetilde{\ell}_{n}^{c}\left(\bm{\Psi}\right)$
$\displaystyle=\sum_{i=1}^{n}\sum_{h=1}^{H}u_{i,h}\left[\log\omega_{h}+\log{f\left(y_{i}\left|\mu_{h}(\bm{x}_{i}),\phi_{h}\right.\right)}\right]$
(3.1)
$\displaystyle\quad-\gamma\sum_{h=1}^{H}{\omega_{h}}{\left\|\bm{\beta}_{h}\right\|_{2}^{2}}-\frac{v}{2}\sum_{h=1}^{H}\sum_{j,k\in
I}\omega_{h}s_{j,k}\lvert{\beta}_{j,h}-{\beta}_{k,h}\rvert.$

Our proposed algorithm consists of two layers of loops to maximize
$\widetilde{\ell}_{n}^{c}\left(\bm{\Psi}\right)$. The outer layer is a classic
EM iteration that repeatedly implements E-steps and M-steps. Assume that
$m=1,2,\dots,M$ represents the the number of EM iterations. The EM algorithm
at the $(m+1)$th iteration is as follows:

(1) E-step. The E step computes the conditional expectation of the function
$\widetilde{\ell}_{n}^{c}\left(\bm{\Psi}\right)$ with respect to $u_{i,h}$,
given the data $(\bm{x}_{i},y_{i})$ and assuming that the current estimate
$\bm{\Psi}^{(m)}$ gives the true parameters of the model. The conditional
expectation (Q-function) is

$\displaystyle Q(\bm{\Psi};\bm{\Psi}^{(m)})$
$\displaystyle=\sum_{i=1}^{n}\sum_{h=1}^{H}\pi_{i,h}^{(m)}\log{\omega_{h}}+\sum_{i=1}^{n}\sum_{h=1}^{H}\pi_{i,h}^{(m)}\log{f\left(y_{i}|\mu_{h}(\bm{x}_{i}),\phi_{h}\right)}$
(3.2)
$\displaystyle-\gamma\sum_{h=1}^{H}{\omega_{h}}{\left\|\bm{\beta}_{h}\right\|_{2}^{2}}-\frac{v}{2}\sum_{h=1}^{H}\sum_{j,k\in
I}\omega_{h}s_{j,k}\lvert{\beta}_{j,h}-{\beta}_{k,h}\rvert,$

where $\pi_{i,h}^{(m)}$ are the expectations of the unobserved $u_{i,h}$ with
the following values:

$\pi_{i,h}^{(m)}=\frac{w_{h}^{(m)}f\left(y_{i}|\mu_{h}^{(m)}(\bm{x}_{i}),\phi_{h}^{(m)}\right)}{\sum\limits_{h^{{}^{\prime}}=1}^{H}w_{h^{{}^{\prime}}}^{(m)}f\left(y_{i}|\mu_{h^{{}^{\prime}}}^{(m)}(\bm{x}_{i}),\phi_{h^{{}^{\prime}}}^{(m)}\right)},$
(3.3)

where
$\mu_{h}^{(m)}(\bm{x}_{i})=g^{-1}(\beta_{0,h}^{(m)}+\bm{x}_{i}\bm{\beta}_{h}^{(m)}).$

(2) M-step. The M step at the $(m+1)$th iteration maximizes the function
$Q(\bm{\Psi};\bm{\Psi}^{(m)})$ with respect to $\bm{\Psi}$. We use a
sequential method to obtain the parameters’ estimates. First, the component
weights $\omega_{h}$ are updated by:

$\omega_{h}^{(m+1)}=\frac{1}{n}\sum_{i=1}^{n}\pi_{i,h}^{(m)},\quad
h=1,2,\dots,H,$ (3.4)

which maximize the leading term of $Q(\bm{\Psi};\bm{\Psi}^{(m)})$ in (3.2).
This is similar to the scheme of Khalili and Chen (2007) and it works well in
our simulations.

Next, we fix $\omega_{h}=\omega_{h}^{(m+1)}$ and maximize
$Q(\bm{\Psi};\bm{\Psi}^{(m)})$ with respect to the other parameters in
$\bm{\Psi}$. Since the pairwise fusion penalty involves a large number of
parameters and is not differentiable with respect to $\beta_{j,h}$, it is
difficult to apply traditional Newton algorithms to realize the M-step
optimization. We further consider re-parameterizing the above optimization
objective function and estimating the parameters within the framework of the
scaled ADMM algorithm444The details for the scaled ADMM can be found in
Appendices A-D.. The inner loops are required in each EM iteration to maximize
the Q-function when using the scaled ADMM. We start from the initial values
$\beta_{0,h}^{(m,0)},\bm{\beta}^{(m,0)}_{h},\phi^{(m,0)}_{h}$ and update the
parameters cyclically for a total of $T$ steps in the $m$th EM iteration. Let
$\bm{z}_{h}=\\{z_{j,k,h}\\}$ be a set of auxiliary variables and
$\bm{r}_{h}=\\{r_{j,k,h}\\}$ be a set of scaled dual variables for the $h$th
component, the scaled ADMM algorithm consists of the following steps:

  * •

Step 1: Initialize $z_{j,k,h}^{(m,0)}=\beta_{j,h}^{(m,0)}$ and
$r_{j,k,h}^{(m,0)}=0$.

  * •

Step 2: At the $(t+1)$th iteration of the scaled ADMM,
$\beta_{0,h}^{(m,0)},\bm{\beta}^{(m,0)}_{h},\phi^{(m,0)}_{h}$ are updated by:

$\displaystyle\mathrm{argmin}_{\beta_{0,h},\bm{\beta}_{h},\phi_{h}}\quad
g(\beta_{0,h},\bm{\beta}_{h},\phi_{h})\triangleq$
$\displaystyle-\sum_{i=1}^{n}\pi_{i,h}^{(m)}\log{f\left(y_{i}|\mu_{h}(\bm{x}_{i}),\phi_{h}\right)}+\omega_{h}\gamma{\left\|\bm{\beta}_{h}\right\|_{2}^{2}}$
$\displaystyle+\frac{\rho}{2}\sum_{j=1}^{p}\sum_{k=1}^{p}\left(z_{j,k,h}^{(m,t)}-\beta_{j,h}+r_{j,k,h}^{(m,t)}\right)^{2}.$
(3.5)

  * •

Step 3: Update the auxiliary variables $\\{z_{j,k,h}^{(m,t+1)}\\}$ by solving
the following optimization problem:

$\mathrm{argmin}_{\bm{z}_{h}}\quad\omega_{h}\frac{v}{2}\sum_{j=1}^{p}\sum_{k=1}^{p}s_{j,k}\left|z_{j,k,h}-z_{k,j,h}\right|+\frac{\rho}{2}\sum_{j=1}^{p}\sum_{k=1}^{p}\left(z_{j,k,h}-\beta_{j,h}^{(m,t+1)}+r_{j,k,h}^{(m,t)}\right)^{2}.$
(3.6)

  * •

Step 4: Update the scaled dual variables $\\{r_{j,k,h}^{(m,t+1)}\\}$, which is
equivalent to the residual between the $\beta_{j,h}^{(m,t+1)}$ and the
auxiliary variable $z_{j,k,h}^{(m,t+1)}$. The estimates can be easily
calculated by:

$r_{j,k,h}^{(m,t+1)}=r_{j,k,h}^{(m,t)}+\left(z_{j,k,h}^{(m,t+1)}-\beta_{j,h}^{(m,t+1)}\right).$
(3.7)

  * •

Step 5: If the termination conditions of the scaled ADMM algorithm are
satisfied, the inner iteration will be stopped in advance, and they are taken
as the estimates of the $(m+1)$th iteration of the EM update. Otherwise, the
scaled ADMM algorithm will go back to Step 2 and continue until reaching the
maximum number of steps $T$.

_Remark_ 1\. We use the standard stopping criterion for the scaled ADMM
algorithm based on primal and dual residuals following Section 3.3 of Boyd et
al. (2011). Specifically, we define the primal and dual residuals
$\delta^{(m,t+1)}_{\text{pri},h}=||\bm{r}_{h}^{(m,t+1)}-\bm{r}_{h}^{(m,t)}||_{2},\
\delta^{(m,t+1)}_{\text{dual},h}=\rho||\bm{z}_{h}^{(m,t+1)}-\bm{z}_{h}^{(m,t)}||_{2}$
for the updates in each component, as well as the criterion
$\delta^{(m,t+1)}_{\text{pri},h}\leq\varepsilon^{\text{pri}}$,
$\delta^{(m,t+1)}_{\text{dual},h}\leq\varepsilon^{\text{dual}}$, where
$\varepsilon^{\text{pri}}$ and $\varepsilon^{\text{dual}}$ are specified as
suggested by Boyd et al. (2011).

(3) Repeat E-step and M-step until the norm difference between the predicted
regression coefficients of the two adjacent EM algorithms
$\lVert\bm{\hat{\mathbf{B}}}^{(m+1)}-\bm{\hat{\mathbf{B}}}^{(m)}\lVert_{F}$ is
less than the given threshold.

The EM-AMDD algorithm is implemented in R programming, as shown in Algorithm
1.

###  3.2 Computational details

Algorithm 1 EM-ADMM algorithm for FMR with covariate clustering.

1:

2:the observed data points
$(y_{1},\bm{x}_{1}),(y_{2},\bm{x}_{2}),\dots,(y_{n},\bm{x}_{n})$, the initial
estimation of the parameters
$\bm{\Psi}^{(0)}=(\omega_{1}^{(0)},\dots,\omega_{H}^{(0)},\beta_{0,1}^{(0)},\dots,\beta_{0,H}^{(0)},\bm{\beta}_{1}^{(0)},\dots,\bm{\beta}_{H}^{(0)},\phi_{1}^{(0)},\dots,\phi_{H}^{(0)})$,
the penalty parameters $\gamma,v,\rho$, the maximum steps $M,T$, the threshold
values
$\varepsilon^{\text{pri}},\varepsilon^{\text{dual}},\varepsilon^{\text{EM}}$,
and the similarity matrix $\mathbf{S}$;

3:the estimated $\bm{\Psi}^{*}$;

4:for $m=0,\dots,(M-1)$ do

5: for $h=1,\dots,H$ do

6: calculate $\pi_{i,h}^{(m+1)}$ by (3.3)-(3.1);

7: update $\hat{\omega}_{h}^{(m+1)}$ by (3.4);

8: initialize $\bm{z}_{h}^{(m,0)}=\\{z_{j,k,h}^{(m,0)}\\}$ and
$\bm{r}_{h}^{(m,0)}=\\{r_{j,k,h}^{(m,0)}\\}$ in the ADMM algorithm;

9: for $t=0,\dots,(T-1)$ do

10: update $\beta_{0,h}^{(m,t+1)},\bm{\beta}_{h}^{(m,t+1)},\phi_{h}^{(m,t+1)}$
by the BFGS algorithm and (A.5)-(A.6);

11: update $\bm{z}_{h}^{(m,t+1)}$ by (A.7)-(A.8);

12: update $\bm{r}_{h}^{(m,t+1)}$ by (3.7);

13: if ($\delta^{(m,t+1)}_{\text{pri}}\leq\varepsilon^{\text{pri}}\ \And\
\delta^{(m,t+1)}_{\text{dual}}\leq\varepsilon^{\text{dual}})\ \text{or}\
t=(T-1)$ then

14: $\hat{\beta}_{0,h}^{(m+1)}\leftarrow\beta_{0,h}^{(m,t+1)}$,
$\hat{\bm{\beta}}_{h}^{(m+1)}\leftarrow\bm{\beta}_{h}^{(m,t+1)}$,
$\hat{\phi}_{h}^{(m+1)}\leftarrow\phi_{h}^{(m,t+1)}$;

15: stop;

16: end if

17: end for ADMM iterations

18: end for all components

19:

20: if
$\lVert\bm{\hat{\mathbf{B}}}^{(m+1)}-\bm{\hat{\mathbf{B}}}^{(m)}\lVert_{F}<\varepsilon^{\text{EM}}\
\text{or}\ m=(M-1)$ then

21: $\bm{\Psi}^{*}\leftarrow\bm{\hat{\Psi}}^{(m+1)}$;

22: stop;

23: end if

24:end for EM iterations

25:return $\bm{\Psi}^{*}$.

(1) Initialization of $\bm{\Psi}^{(0)}$. Although numerical methods like EM
and ADMM have good convergence properties for convex optimization problems,
the choice of initial values remains important, particularly for problems with
a large number of parameters to estimate. Using an uninformative prior can
lead to suboptimal results, affecting the algorithm’s convergence speed and
stability. Here, we propose a method that combines sample clustering with
regularized GLM regression to address this issue. The identification of
subgroups of the outcome corresponds to the initial component weights of FMR
models. Specifically, we perform K-means clustering on observations $y_{i}$ to
derive the prior weights of the components. This approach is based on the idea
that mixing is equivalent to classifying the latent risk level of each
individual. Using the clustering result, we partition the sample into $H$
subgroups and train a Ridge GLM model on each subgroup to obtain initial
estimates of $\beta_{0,h}$, $\bm{\beta}_{h}$, and $\phi_{h}$, respectively.
This can provide a starting point for the subsequent EM or ADMM optimization.

(2) Selection of the tuning parameters. The three tuning parameters in our
model, $\gamma$, $v$, and $\rho$, each have a distinct role. Specifically,
$\gamma$ controls the extent of shrinkage in the regression coefficients, $v$
regulates the penalty on the distance between similar covariates, and $\rho$
affects the convergence efficiency of the ADMM algorithm. For our study, we
set $\rho=1$ and determined the optimal values of $\gamma$ and $v$ through
cross-validation on the training set. As we are primarily interested in the
impact of the $\ell_{2}$ norm penalty and covariate clustering penalty, we
evaluated the tuning parameter selection criteria using log-likelihood value
or mean square error (MSE) criterion, given that our response variable of
interest is continuous healthcare costs. To provide insight into the
computational requirements of our proposed EM-ADMM algorithm, we implemented
the simulation setup outlined in Section 4, where $n=2000$ and $p=10$. On a
computer equipped with a 3.61 GHz Intel Core-I7 CPU and 32 GB RAM, each EM-
ADMM iteration took an average of 1.3 seconds to run. Matrix operations can be
utilized to expedite the algorithmic speed.

(3) Identifying the covariates’ clusters. Although the objective function
(2.4) theoretically guarantees that the difference between the coefficients of
similar covariates will converge to zero, this result is difficult to achieve
in practice due to limitations of numerical algorithms. Therefore, we follow
the strategy proposed in Chi and Lange (2015), which places an edge between
covariates $j$ and $k$ if $z_{j,k,h}$ is equal to $z_{k,j,h}$. Then, two
covariates sharing the same edge can be classified into the same cluster in
the $h$th component. This condition is established when $\theta=0.5$ in (A.7).

(4) Number of iterations and stopping criterion. We set $M=10$ and $T=100$,
which means that the EM algorithm is executed at most 10 times and each M-step
involves a maximum of 100 ADMM iterations. When determining the early
termination conditions of the ADMM and EM algorithms, it is essential to
consider the magnitude and dimension of the regression coefficients. In cases
where the dimensionality of the covariates is high, we recommend setting the
ADMM thresholds to $\varepsilon^{\text{pri}}=0.05$ and
$\varepsilon^{\text{dual}}=0.05$, and the EM threshold to
$\varepsilon^{\text{EM}}=0.01$. We have found that these threshold values
strike a balance between computational speed and convergence. This balance is
achieved by monitoring the value of the log-likelihood function in subsequent
numerical simulations and empirical studies.

##  4 Simulation

A simulation study will be conducted to verify the effect of the proposed
model on clustering and coefficient estimation for multiple groups of
covariates. The observations are assumed to be generated by a mixture
distribution consisting of two components, where each component is a Gamma
distribution with a mean parameter that depends on the covariates:

$Y_{i}\sim\omega_{1}\text{Gamma}(y_{i}\left|\mu_{1}(\bm{x}_{i}),\phi_{1}\right.)+\omega_{2}\text{Gamma}(y_{i}\left|\mu_{2}(\bm{x}_{i}),\phi_{2}\right.),\quad
i=1,\dots,n,$

where $\mu_{h}(\bm{x}_{i})$ and $\phi_{h}(h=1,2)$ respectively denote the mean
and dispersion parameters, and $\omega_{1}+\omega_{2}=1$. The mean parameter
$\mu_{h}(\bm{x}_{i})$ can be expressed as a function related to the covariates
$\bm{x}_{i}$ through a log link function, specifically,
$\log\mu_{h}(\bm{x}_{i})=\beta_{0,h}+\bm{x}_{i}^{T}\bm{\beta}_{h}$.

The true values of the weights, regression coefficients and dispersion of the
mixture Gamma model are set to be as follows:

$\displaystyle\text{Component}\ 1$ $\displaystyle:\omega_{1}=0.7,\
\beta_{0,1}=1,\
\bm{\beta}_{1}=(\overbrace{\beta_{1,1},...,\beta_{5,1}}^{-0.1},\overbrace{\beta_{6,1},...,\beta_{10,1}}^{-0.2}),\
\phi_{1}=0.2;$ $\displaystyle\text{Component}\ 2$
$\displaystyle:\omega_{2}=0.3,\ \beta_{0,2}=2,\
\bm{\beta}_{2}=(\overbrace{\beta_{1,2},...,\beta_{5,2}}^{0.1},\overbrace{\beta_{6,2},...,\beta_{10,2}}^{0.2}),\
\phi_{2}=0.1;$

which is equivalent to dividing the 10 covariates of the $i$th observation
$\bm{x}_{i}=(1,\overbrace{x_{i1},...,x_{i5}}^{\text{cluster
1}},\overbrace{x_{i6},...,x_{i10}}^{\text{cluster 2}})$ into 2 clusters. Each
covariate in the same cluster shares a common effect on the outcome data. For
example, in component 1, the coefficient of covariate $k$ is 0.1 for
$k=1,\dots,5$, and 0.2 for $k=6,\dots,10$. The covariates have completely
opposite regression coefficients in the two components, reflecting the
challenges caused by unobserved heterogeneous risk factors in healthcare cost
modeling.

We generate the dataset with a sample size $n=1000$ during each run of the
simulation. The covariates $\bm{x}_{i}$ are assumed to be generated by a
multivariate normal distribution $\text{MN}(\mathbf{0},\mathbf{\Sigma})$ with
the mean of $\mathbf{0}$ and the covariance matrix of $\mathbf{\Sigma}$.
$\mathbf{\Sigma}$ is a $p\times p$ block diagonal matrix, with the elements as
follows:

$\Sigma_{kk^{\prime}}=\begin{cases}0.04&\text{if}\quad k=k^{\prime},\\\
0.04\varrho&\text{if}\quad 1\leq k\leq 5,1\leq k^{\prime}\leq 5,k\neq
k^{\prime},\\\ 0.04\varrho&\text{if}\quad 6\leq k\leq 10,6\leq k^{\prime}\leq
10,k\neq k^{\prime},\\\ 0&\text{otherwise},\end{cases}$

where $\varrho$ denotes the correlation coefficient. We can vary the value of
$\varrho$ from the grid $\\{0.1,0.5,0.9\\}$ to explore the influence of
covariate similarity on our model. Thus, the observations $Y_{i}$ for
$i=1,\dots,n$ are generated according to the following scenario: first, a
random number $u_{i}$ is drawn from a uniform distribution over the interval
$(0,1)$. If $u_{i}$ is less than $\omega_{1}$, $Y_{i}$ is drawn from the first
component, specifically
$\text{Gamma}(y_{i}\left|\mu_{1}(\bm{x}_{i}),\phi_{1}\right.)$; otherwise, it
is drawn from the second component, specifically
$\text{Gamma}(y_{i}\left|\mu_{2}(\bm{x}_{i}),\phi_{2}\right.)$.

For each simulated dataset, we fit the proposed model with a wide range of
tuning parameter values $v$555The other two tuning parameters are fixed at
$\gamma=0.001$ and $\rho=1$.. Figure 1 displays the solution paths against $v$
for clustering covariates under the three simulated samples of
$\varrho=0.1,0.5,0.9$. It is evident that if the effects of the covariates on
the response variable exhibit different clustering effects in different
subsamples, the proposed model will correctly cluster the covariates by
identifying the same corresponding regression coefficients in each component.
Our model also tends to cluster covariates when they have a certain
correlation (e.g., $\varrho=0.5$), and it shows better identification of the
main cluster with greater component weight in FMR. In Table 1, we report the
correctly clustered proportion (CCP) of the covariates among all observations
of each latent component, which measures the agreement between the true and
estimated clusters, and a value close to 1 is preferred. As expected, CCP-1
and CCP-2 both increase as the tuning parameter $v$ increases, especially when
the true correlation of the covariates is high. The impact of covariate
clustering is not apparent when the correlation is weak.

Figure 1: The solution paths for the covariates clustering against $v$ based
on random samples with various values of $\varrho$. Left panel and Right panel
display the estimated coefficients in component 1 and 2 respectively. The
values of $\varrho=0.1,0.5,0.9$ are presented in the Top, Middle, and Bottom
panels, respectively. The true valuesof regression coefficients are marked as
the black solid lines.  Table 1: The correct clustering proportion (CCP) in
each component for Gamma mixture regressions in the case that the covariates
are assumed to be grouped into two clusters.

$v$ | $\varrho=0.1$ | $\varrho=0.5$ | $\varrho=0.9$
---|---|---|---
CCP-1 | CCP-2 | CCP-1 | CCP-2 | CCP-1 | CCP-2
0 | 20.00% | 20.00% | 20.00% | 20.00% | 20.00% | 20.00%
4 | 40.00% | 36.00% | 64.00% | 36.00% | 52.00% | 44.00%
8 | 44.00% | 36.00% | 88.00% | 56.00% | 88.00% | 92.00%
12 | 44.00% | 36.00% | 92.00% | 76.00% | 100.00% | 100.00%
16 | 44.00% | 36.00% | 96.00% | 92.00% | 100.00% | 100.00%
20 | 44.00% | 40.00% | 100.00% | 100.00% | 100.00% | 100.00%

  * •

* CCP-1 and CCP-2 denote the correct clustering proportion in component 1 and 2 respectively.

We further conduct an investigation into the performance of the EM-ADMM
algorithm in estimating parameters with respect to sample size $n$ and tuning
parameter $v$, assuming high correlation of covariates in groups
($\varrho=0.9$). The results of this evaluation are presented in Table 2,
where bias and MSE were computed based on 100 repetitions for each of the
sample sizes ($n=100,500,1000$). Bias is defined as the difference between the
average estimate and the true value of the regression coefficient, with a
value close to zero indicating good performance. MSE measures the average
squared distance from the estimate to the true value, which is a commonly used
metric for assessing prediction performance. Table 2 indicates that the model
provides reasonable estimates for each coefficient, as long as a sufficient
penalty parameter $v$ is added, regardless of the small or large sample size.
The accuracy of the parameter estimation results and the stability of the
optimization algorithm are also reflected in the table. As the amount of
training data increases, the estimation of the model becomes more accurate,
and the magnitude of the bias and the MSE becomes almost negligible in
comparison with the true values of the regression coefficients.

Table 2: Bias and MSE for estimated regression coefficients for Gamma mixture models in the case of high correlation level of $\varrho=0.9$. |  | $n=100$ | $n=500$ | $n=1000$
---|---|---|---|---
|  | $v=0$ | $v=20$ | $v=0$ | $v=20$ | $v=0$ | $v=20$
|  | Bias | MSE | Bias | MSE | Bias | MSE | Bias | MSE | Bias | MSE | Bias | MSE
Com1 | $\beta_{1,1}$ | -0.031 | 0.045 | -0.003 | 0.000 | -0.005 | 0.005 | 0.001 | 0.000 | 0.001 | 0.004 | -0.001 | 0.000
$\beta_{2,1}$ | -0.021 | 0.047 | -0.002 | 0.000 | -0.006 | 0.006 | 0.000 | 0.000 | -0.003 | 0.003 | -0.001 | 0.000
$\beta_{3,1}$ | 0.023 | 0.055 | -0.001 | 0.000 | 0.010 | 0.007 | 0.000 | 0.000 | -0.004 | 0.003 | 0.000 | 0.000
$\beta_{4,1}$ | 0.013 | 0.051 | -0.002 | 0.000 | -0.004 | 0.007 | -0.001 | 0.000 | -0.010 | 0.003 | 0.000 | 0.000
$\beta_{5,1}$ | 0.006 | 0.033 | -0.003 | 0.000 | 0.003 | 0.005 | -0.001 | 0.000 | 0.013 | 0.003 | 0.000 | 0.000
$\beta_{6,1}$ | -0.001 | 0.046 | 0.000 | 0.000 | -0.003 | 0.008 | -0.002 | 0.000 | -0.007 | 0.003 | -0.001 | 0.000
$\beta_{7,1}$ | -0.021 | 0.041 | 0.000 | 0.000 | -0.016 | 0.005 | 0.000 | 0.000 | 0.009 | 0.003 | 0.002 | 0.000
$\beta_{8,1}$ | 0.001 | 0.038 | 0.000 | 0.000 | 0.007 | 0.005 | -0.002 | 0.000 | 0.002 | 0.003 | -0.001 | 0.000
$\beta_{9,1}$ | 0.012 | 0.037 | -0.001 | 0.000 | -0.001 | 0.005 | -0.001 | 0.000 | -0.009 | 0.003 | -0.002 | 0.000
$\beta_{10,1}$ | -0.001 | 0.036 | 0.002 | 0.000 | 0.008 | 0.005 | 0.000 | 0.000 | 0.001 | 0.003 | -0.001 | 0.000
Com2 | $\beta_{1,2}$ | -0.082 | 0.535 | -0.003 | 0.002 | 0.022 | 0.078 | -0.006 | 0.000 | 0.012 | 0.028 | 0.000 | 0.000
$\beta_{2,2}$ | -0.044 | 0.344 | -0.002 | 0.002 | -0.034 | 0.053 | -0.005 | 0.000 | 0.028 | 0.030 | 0.001 | 0.000
$\beta_{3,2}$ | 0.060 | 0.475 | -0.003 | 0.002 | 0.019 | 0.044 | -0.005 | 0.000 | -0.025 | 0.030 | 0.000 | 0.000
$\beta_{4,2}$ | 0.053 | 0.427 | -0.003 | 0.002 | -0.028 | 0.051 | -0.006 | 0.000 | -0.012 | 0.029 | -0.001 | 0.000
$\beta_{5,2}$ | -0.031 | 0.375 | 0.000 | 0.002 | -0.012 | 0.064 | -0.005 | 0.000 | -0.011 | 0.027 | -0.002 | 0.000
$\beta_{6,2}$ | 0.024 | 0.396 | -0.011 | 0.003 | -0.006 | 0.052 | -0.004 | 0.000 | 0.007 | 0.023 | 0.001 | 0.000
$\beta_{7,2}$ | -0.007 | 0.445 | -0.014 | 0.003 | -0.007 | 0.048 | -0.004 | 0.000 | 0.004 | 0.026 | 0.001 | 0.000
$\beta_{8,2}$ | -0.048 | 0.415 | -0.012 | 0.003 | -0.007 | 0.059 | -0.003 | 0.000 | -0.004 | 0.026 | -0.001 | 0.000
$\beta_{9,2}$ | -0.099 | 0.505 | -0.012 | 0.002 | -0.003 | 0.052 | -0.003 | 0.000 | -0.025 | 0.027 | 0.000 | 0.000
$\beta_{10,2}$ | 0.053 | 0.375 | -0.012 | 0.003 | 0.002 | 0.053 | -0.003 | 0.000 | 0.014 | 0.027 | 0.001 | 0.000

##  5 Empirical analysis

###  5.1 Experimental setup

In this section, we present two practical examples to demonstrate the proposed
covariate clustering regression model. The first example uses a dataset on
diabetes expenses from several hospitals in China, which includes patient
characteristics and questionnaire data. The aim is to investigate the
relationship between annual expenses and individuals’ knowledge of diabetes.
The second example uses a short-term medical claim dataset collected from a
Chinese insurance company. The objective is to explore the clustered effect of
the ICD (International Classification of Diseases) categories in the subgroups
of the medical insurance claim data. To prevent model overfitting, we divide
the dataset $\mathcal{L}(y_{i},\bm{x}_{i}){0\leq i\leq n}$ into a training
set, which is used for model fitting and hyperparameter tuning, and a testing
set, which is used for out-of-sample analysis only.

A comparative study is conducted for model comparison in terms of goodness of
fit, predictive accuracy, and discriminative ability of samples of different
risk in the testing dataset. Two commonly used metrics to evaluate the
goodness of fit of a model are the negative log-likelihood (NLL) and the
pseudo R2. The Gamma distribution is chosen as the underlying distribution
assumption for all models, so the density function $f(y)$ in (2.3) has the
form

$f\left(y\left|\mu,\phi\right.\right)=\frac{y^{\frac{1}{\phi^{2}}-1}\exp\left(-\frac{y}{\mu\phi^{2}}\right)}{(\mu\phi^{2})^{\frac{1}{\phi^{2}}}\
\Gamma({1}/{\phi^{2}})},$ (5.1)

where $\Gamma(\alpha)=\int_{0}^{\infty}t^{\alpha-1}e^{-t}dt$ is the gamma
function, and NLL can be defined as $-\ell_{n}(\bm{\Psi})$ with a set of
estimated parameters. The pseudo R2 measures the gap of log-likelihood between
the saturated and predictive model. Specifically, given the sample
observations $\bm{y}=({y}_{1},\dots,{y}_{n})$ and the predictions
$\hat{\bm{y}}=(\hat{y}_{1},\dots,\hat{y}_{n})$, where
$\hat{y}_{i}=\sum_{h=1}^{H}\hat{\omega}_{h}\hat{\mu}_{h}(\bm{x}_{i})$, the
pseudo R2 under the Gamma distribution assumption can be expressed as:

$\text{R}^{2}=1-\frac{\sum_{i=1}^{n}\left[\log(y_{i}/\hat{y}_{i})-(y_{i}-\hat{y}_{i})/\hat{y}_{i}\right]}{\sum_{i=1}^{n}\log(y_{i}/\bar{y})},$
(5.2)

where $\bar{y}$ is the sample mean of $\bm{y}$. A low value of NLL and a high
value of R2 are preferred for the model. Regarding predictive accuracy, we
employ two metrics: mean squared error (MSE) to evaluate the accuracy of point
estimation, and mean continuous ranked probability score (MCRPS) to evaluate
the accuracy of distribution estimation. While MSE is easy to calculate, CRPS
measures how well the overall predicted distribution takes values around the
true observation. Its analytical form is given by Gneiting and Raftery (2007)
as:

$\text{CRPS}(y_{i},\hat{F})=\int_{\mathbb{R}}\left[\hat{F}(y_{i}|\bm{x}_{i})-\text{I}(y_{i}\leq
z)\right]^{2}dz,$ (5.3)

where $\hat{F}(\cdot)$ is the estimated cumulative distribution function
corresponding to (2.1), and $\text{I}(\cdot)$ is an indicative function. A
predicted distribution with a smaller CRPS would logically be more
concentrated around the true observation. We then average the CRPS values of
all samples to obtain the model’s performance on the entire dataset. Finally,
we use the two-way lift metric to assess the model’s ability to discriminate
risk across samples. To calculate the lift, we bin the samples into ten equal
groups according to their predicted losses, from low to high. We then divide
the average response of the top-decile group by that of the bottom-decile
group (Lee, 2021). A higher lift means that the model can better distinguish
between high-risk and low-risk samples.

###  5.2 Diabetes expense dataset

The first dataset comprises healthcare costs data for 5,961 diabetes patients
collected from several large hospitals in China during the year 2018. Each
observation records the positive annual expense of a patient that covers
treatment and hospitalization costs incurred due to diabetes. The dataset
contains patient-specific covariates related to the healthcare costs,
primarily demographic characteristics such as gender, age, monthly income,
payment method, and basic health status such as smoking, drinking, duration of
diabetes, body mass index, checkup values of blood pressure and blood glucose.
The dataset also includes information on complications of diabetes, such as
chronic complications (including hypertension and vessel disease) and acute
complications, as well as medication use, including whether oral hypoglycemic
drugs or insulin are used. These variables are aggregated to form 28
traditional risk factors. Additionally, a questionnaire survey is conducted
for each patient to evaluate their knowledge of diabetes, which consists of 26
main multiple-choice questions, with a total of 111 options. The questionnaire
survey records whether each option is answered correctly, denoted as binary
variables of 0 and 1.

We aim to investigate the association between the individuals’ knowledge of
diabetes and diabetes healthcare cost. However, directly modeling the original
111 covariates is not advisable, as the information reflected by different
options of each question, or even different questions, may overlap to some
extent. Consequently, we will only cluster the 111 covariates of the
questionnaire in the follow-up experiments and retain the traditional
covariates without further processing. Synchronized covariate clustering in
regression modeling is a viable technique for reducing dimensionality and
mining the textual information of the questionnaire. We adopt a data-driven
approach to construct the similarity matrix $\mathbf{S}={s_{j,k}}_{p\times p}$
by computing the cosine similarity of covariate $j$ and covariate $k$ in a
training set, while more sophisticated methods of defining covariate
similarity are worth exploring in future research.

Table 3 provides an overview of the 139 covariates detailing the risk factors.
Figure 2 illustrates the empirical distribution of diabetes expenses. The log-
log plot indicates that the largest 10% of the expense tail exhibits a linear
trend, implying that the response variable is heavy-tailed. Additionally, the
density function is highly skewed and, to some extent, multi-modal; however,
the subsequent peaks are not significant.

Figure 2: Diabetes expense dataset: Left panel: histogram of the patients’ healthcare costs; Middle panel: empirical density (upper-truncated at 75,000); Right panel: log-log plot of observed expenses. Table 3: Diabetes expense dataset: description of variables. Variable | Type | Description
---|---|---
healthcare costs | Continuous | Annual healthcare costs of an inpatient for diabetes: 3-200,000 CNY.
Gender | Categorical | Binary: male or female.
Age | Continuous | Integer number: 16-94 years.
Education level | Categorical | Six types: from below elementary school to undergraduate and above.
Marital status | Categorical | Binary: married or not.
Occupation | Categorical | Four types: on-the-job, retired, laid-off, other.
Income | Continuous | Monthly income of the patient: 24-100,000 CNY.
Pay way | Categorical |  | Payment way of healthcare costs, four types: self-pay, public-pay,
---
medical insurance, other.
Smoking | Categorical | Binary: smoking or not.
Drinking | Categorical | Binary: drinking or not.
Duration of diabetes | Continuous | How long since diabetes was diagnosed: 1-84.7 years.
Physical indicators | Continuous |  | A total of 6 covariates, reflecting the basic health status of the inpatient.
---
These covariates include: body mass index (12.80-60.44), systolic blood
pressure (60-240mmHg), diastolic blood pressure (30-180mmHg), fasting
plasma glucose (0.56-43.20mmol/L), 2-hour postprandial glucose (0.9-55
mmol/L), and hemoglobin A1c (3-18mmol/L). All are continuous.
Chronic complications | Categorical |  | A total of 7 covariates, recording whether the inpatient has chronic
---
complications of diabetes (i.e., hypertension, vessel disease, diabetic
nephropathy, diabetic neuropathy, diabetic eye disease, diabetic foot
disease, and others). All are binary variables, with “1” being diseased.
Acute complications | Categorical |  | A total of 3 covariates, recording whether the inpatient has acute comp-
---
lications of diabetes (i.e., diabetic ketoacidosis, hypoglycemia, and lactic
acidosis). All are binary variables, with “1” being diseased.
Oral medication | Categorical | Binary: whether oral hypoglycemic drugs are used or not.
Insulin medication | Categorical | Binary: whether insulin is used or not.
Questionnaire data | Categorical |  | A total of 111 covariates, measuring the individual’s understanding of
---
diabetes-related knowledge. All are binary variables: a certain covariate
is recored as “1” if the corresponding option is answered correctly.

Table 4 shows the out-of-sample performance of competing models, covering a
comprehensive comparison among conventional GLMs, ridge GLMs and covariate
clustering GLMs, as well as the corresponding mixture models. Specifically, M1
to M4 represent the conventional GLM (without $\ell_{2}$-penalty and
similarity penalty), the ridge GLM (without similarity penalty), the mixture
GLM (without $\ell_{2}$-penalty and similarity penalty), and the mixture ridge
GLM (without similarity penalty), respectively. M5-M8 are similar to M1-M4,
but the last four models all incorporate our covariate clustering framework
(i.e., with similarity penalty). To speed up the convergence of the algorithm,
the response variable is transformed into $y^{*}=y/10,000$, all continuous
covariates are normalized, and the other (multiple) categorical covariates are
treated with dummy variables. The dataset is randomly split into the training
and testing sets with a proportion of 80:20, leading to 4,768 observations for
training and 1,193 observations for testing. Regarding tuning parameters
selection, when the $\ell_{2}$ norm regularization term is considered,
$\gamma$ is tuned from the three values ${0.5,1.0,5.0}$; otherwise,
$\gamma=0.01$ is fixed to ensure the strict convexity of the objective
function. When the regularization term for covariate clustering is added,
$\rho=1.0$ is fixed, and $v$ is tuned from the grid values
$\\{0.5,1.0,\dots,4.5,5.0\\}$; otherwise, $\rho=0.001$, and $v=0$ are fixed.
We select the tuning parameters that can achieve the minimum NLL on the
validated set for each model through further subdividing a validated part from
the training set.

Table 4: Diabetes expense dataset: out-of-sample performance of the eight
competing models.

Models | Hyperparameters | Out-of-sample Performance
---|---|---
H | $\gamma$ | $\rho$ | $v$ | NLL | Pseudo R2 | MSE(10,000) | MCRPS | Lift
M1 | 1 | 0.01 | 0.001 | 0 | 733.718 | 15.83% | 1.384 | 0.363 | 3.588
M2 | 1 | 0.5 | 0.001 | 0 | 733.743 | 15.83% | 1.384 | 0.363 | 3.588
M3 | 2 | 0.01 | 0.001 | 0 | 801.783 | 11.97% | 1.373 | 0.356 | 3.716
M4 | 2 | 5.0 | 0.001 | 0 | 812.429 | 11.08% | 1.383 | 0.357 | 3.606
M5 | 1 | 0.01 | 1.0 | 0.5 | 716.809 | 17.51% | 1.388 | 0.359 | 4.250
M6 | 1 | 1.0 | 1.0 | 0.5 | 717.408 | 17.44% | 1.389 | 0.359 | 4.210
M7 | 2 | 0.01 | 1.0 | 3.5 | 780.983 | 14.04% | 1.368 | 0.354 | 3.723
M8 | 2 | 1.0 | 1.0 | 2.5 | 793.535 | 12.29% | 1.376 | 0.354 | 4.064

  * •

*The optimal values of metrics are in bold.

Based on the results presented in Table 4, several main conclusions can be
drawn. Firstly, the application of the covariate clustering method
significantly improves the predictive performance of diabetes costs in this
small dataset characterized by a large number of highly correlated risk
factors. The comparison between models M5 to M8 and models M1 to M4 reveals
that most of the metrics of the former set are superior to those of the
latter, with the exception of the MSEs of M5 and M6, which are slightly worse
than those of M1 and M2. Notably, the covariate clustering model not only
accurately estimates the expected value (point estimates) but also effectively
characterizes the overall distribution of the target outcome, as evidenced by
the improvement in goodness of fit (i.e., NLL and pseudo R2) and the ability
to discriminate extreme observations (i.e., lift). This phenomenon is
consistent with our expectations that the covariate clustering method is
equivalent to aggregating the effects of multiple related risk factors, which
eliminates the random disturbance and mutual interference among the
covariates, ultimately enhancing the data.

Secondly, once the covariates with similar or repeated effects have been
clustered during the regression process, the original $\ell_{2}$ norm penalty
for the compression and selection of covariates is no longer necessary. The
results indicate that the performance of models M2, M4, M6, M8 is not as
optimal as that of models M1, M3, M5, M7. As previously mentioned,
dimensionality reduction can also be achieved by clustering similar
covariates, rendering the use of the $\ell_{2}$ norm penalty unnecessary in
subsequent experiments.

Thirdly, when comparing mixture and non-mixture models, model M5 exhibits the
best performance in terms of goodness of fit and risk discrimination ability,
while M7 performs best in terms of prediction accuracy. Although non-mixture
models can effectively capture the characteristics of diabetes expenses data
where the heterogeneity is insignificant, mixture models can generally capture
the majority values of the sample and have advantages in point estimation.

As demonstrated previously, our proposed framework effectively integrates
prior similarity with model performance to generate meaningful posterior
clusters and regression coefficients for the covariates. We employ network
analysis to examine the relationships between the covariates in the resulting
model. Figure 3 depicts the unweighted and undirected graphs. It is worth
noting that the distance between non-connected clusters does not reflect their
similarity. In component 1, the covariates are predominantly grouped into one
main cluster (located in the middle) and two minor clusters (in the upper-
right and lower-right regions), while the remaining covariates exhibit unique
effects that cannot be aggregated. Component 2, however, yields more distinct
clusters, allowing for better differentiation of the covariates’ effects. This
finding is in line with the prior similarity analysis. In other words, the
model clustering most covariates in subgroup 1 (as in component 1, which
accounts for approximately $\hat{\omega}_{1}=0.96$ of the mixture model’s
weight), while it fine-tunes the clustering effects on subgroup 2 (as in
component 2, $\hat{\omega}_{2}=0.04$).

Furthermore, the utilization of colored backgrounds provides a more nuanced
representation of clustering outcomes, where community formation is determined
by the edge betweenness of each vertex (Girvan and Newman, 2002). Moreover, we
have identified and compiled a list of significant covariates for large
communities, enabling further exploration of their internal relationships. For
instance, component 1’s largest cluster can be subdivided into four
communities, with one of them containing “V4” ,“V5” ,“V77”, and “V79,” which
relate to fundamental knowledge and the detrimental effects of smoking on
diabetes patients. The exhaustive description of these variables is available
in Appendices D.

(a) covariate clustering results in component 1  (b) covariate clustering
results in component 2

Figure 3: Diabetes expense dataset: clustering results of the proposed model:
the vertex is a covariate, the edge is an association between two covariates,
and the covariates can be grouped into a single cluster if a connected path
exists between them; the colored background depicts a covariate’s “community”
with finer granularity than the “cluster”.

The comprehensive findings underscore the limitations of implementing
unsupervised clustering and the conventional regression methods in isolation.
Our framework concurrently addresses both tasks to attain more precise results
for both regression and covariate clustering. While we have solely conducted
elementary analyses based on the covariates’ network structure, adopting more
semantic methods in the future can enhance the similarity matrix’s definition
and yield more comprehensive outcomes regarding patients’ diabetes self-
awareness.

###  5.3 Short-term medical claim dataset

As a second example, we examine a medical insurance dataset provided by a
major insurance company in China, which includes individual inpatient claim
data for a short-term medical insurance plan covering the period from 2014 to
2016. The dataset comprises 19,110 policies, and each policy includes the
positive inpatient claim amount as well as several covariates such as Gender,
Age, SSCoverage, HospitalDays, ClaimType, and International Classification of
Diseases (ICD) code. Table 5 displays the summary statistics of these
variables.

Figure 5 shows how the covariates are distributed in the medical insurance
dataset. We model the claim amounts by incorporating covariates through the
two-component mixture Gamma regression models for covariates clustering,
illustrating how clustering the ICD dummy variables can improve the medical
claims prediction. For our analysis, we split the dataset into a training set
and a testing set in proportion 70:30. To increase the convergence of the
algorithm, two continuous covariates (Age and HospitalDays) are standardized.
The ICD variables are transformed into 20 dummy variables. We consider the
non-informative similar matrix for these ICD dummy variables by setting
$S_{kk^{\prime}}=0.5$ for $k,k^{\prime}\in\left\\{1,\dots,20\right\\}$. For
non-similarity covariate $i\neq k$ and $j\neq k^{\prime}$, we set $S_{ij}=0$.
The regression coefficients are estimated by using proposed EM-ADMM with
$\rho=1$ and $\gamma=0.01$. The optimal tuning parameter for covariate
clustering is $v=4.5$, which is selected from
$v\in\left\\{0,0.5,1,\dots,5\right\\}$ based on the minimum NLL criteria on
training data.

Figure 4 shows the estimated coefficients for two components respectively in
mixture regression model with different values of $v$. An ideal covariate
clustering method should achieve high prediction accuracy with few clusters,
as fewer clusters lead to more compact models and potentially easier
interpretability. We find that the number of clustering groups for the ICD
covariates in the first second component tends to decrease as the tuning
parameter $v$ increases, indicating that the covariate clustering method for
mixture regression leads to significantly more compact models than the mixture
model without covariate clustering, especially in the first component.
Specifically, the proposed method makes it possible to classify ICD18 and ICD4
into one category by taking into account the outcome’s information. However,
since the weight parameter in the second component is estimated for the range
of $w_{2}\in\left[0,0.05\right]$, which gives less weight to the penalty of
the log-likelihood function in the second component, the proposed method has
different performance in each subject subgroup.

Figure 4: Medical claim dataset: the estimated regression coefficients for two components in mixture Gamma regression model when $v$ varies from 0 to 5 in the training dataset. The optimal $v=4.5$ is marked as the vertical line. The estimated weights for the two components are $\hat{w}_{1}=95.33\%$ and $\hat{w}_{2}=4.67\%$ respectively.  Table 5: Medical claim dataset: description of variables. Variables | Type | Description
---|---|---
Claim amount | Continuous | inpatient (positive) claim amount of a patient ¥3-200,000
Age | Continuous | inpatient’s age: 18-67
Gender | Categorical | gender of a inpatient (2 levels): male, female
SSCoverage | Categorical | social medical insurance or public medical insurance coverage (2 levels): yes or no
HospitalDays | Continuous | the duration of hospitalization of a inpatient: 0-184 days
ClaimType | Categorical |  | reasons for claiming the medical insurance claim (3 levels):
---
medical treatment from disease (MTD),
medical treatment from accident (MTA), other
ICD | Categorical |  | ICD-10 (20 levels):
---
1: A00-B99, Certain infectious and parasitic diseases
2: C00-D49, Neoplasms
3: D50-D89, Diseases of the blood and blood-forming organs and certain
disorders
involving the immune mechanism
4: E00-E89, Endocrine, nutritional and metabolic diseases
5: F01-F99, Mental, Behavioral and Neurodevelopmental disorders
6: G00-G99, Diseases of the nervous system
7: H00-H59, Diseases of the eye and adnexa
8: H60-H95, Diseases of the ear and mastoid process
9: I00-I99, Diseases of the circulatory system
10: J00-J99, Diseases of the respiratory system
11: K00-K95, Diseases of the digestive system
12: L00-L99, Diseases of the skin and subcutaneous tissue
13: M00-M99, Diseases of the musculoskeletal system and connective tissue
14: N00-N99, Diseases of the genitourinary system
15: O00-O9A, Pregnancy, childbirth and the puerperium
17: Q00-Q99, Congenital malformations, deformations and chromosomal
abnormalities
18: R00-R99, Symptoms, signs and abnormal clinical and laboratory findings,
not elsewhere classified
19: S00-T88, Injury, poisoning and certain other consequences of external
causes
20: V00-Y99, External causes of morbidity
21: Z00-Z99, Factors influencing health status and contact with health
services
Figure 5: Medical claims case: relative frequency of the covariates Gender,
Age, SSCoverage, ICD, HospitalDays and ClaimType.

To illustrate the superiority of the proposed models, we consider fix
competing regression models: the GLMs without ICD (M1), the GLMs with ICD
dummies (M2), the GLMs with ICD clustering (M3), the mixture model without ICD
(M4), the mixture model with ICD dummies (M5), the mixture model with ICD
clustering (M6). Table 6 reports the NLL, Pseudo $R^{2}$, MSE, MCRPS, and Lift
metrics of the different fitted models on the testing data. We select the
optimal tuning parameter $v=0.1$ of M3 from
$v\in\left\\{0,0.1,0.2,\dots,5\right\\}$ based on the minimum NLL criteria on
the training data. Our results show that the two-component mixture models
(M4-M6) outperform the traditional GLMs (M1-M3) as evidenced by their lower
NLL values. This suggests the presence of unobserved latent factors causing
heterogeneity in claim amounts. The QQ-plots of the quantile residuals
provided in Figure 6 support this finding, indicating a preference for the
two-component mixture models in capturing the multi-modal and heavy-tailed
features of medical claim amounts. Furthermore, the covariate clustering
method for mixture models improves model performance in terms of out-of-sample
MSE, with M6 having the lowest MSE. However, as expected, M6 does not perform
well in terms of Pseudo $R^{2}$ and MCRPS metrics, while showing better
ability to distinguish good risks from bad risks in terms of the lift metric.
Figure 7 illustrates this finding with lift plots for selected candidates.
Overall, our results demonstrate the effectiveness of the proposed models in
capturing the heterogeneity in medical claim amounts and improving model
performance.

Figure 6: Medical claim dataset: QQ-plot for quantile residuals on the testing
data. Table 6: Medical claim dataset: out-of-sample performance for the six
competing models.

Models | Out-of-sample Performance
---|---
NLL | Pseudo R2 | MSE ($10^{11}$) | CRPS | Lift
M1 | 54420.74 | 0.264 | 1.15 | 0.0067427 | 5.89
M2 | 54426.70 | 0.262 | 1.10 | 0.0067428 | 5.74
M3 | 54426.21 | 0.262 | 1.15 | 0.0067428 | 5.75
M4 | 53935.07 | 0.260 | 2.37 | 0.0067474 | 5.83
M5 | 53945.57 | 0.259 | 2.88 | 0.0067475 | 5.91
M6 | 54007.50 | 0.258 | 0.18 | 0.0067465 | 6.16

  * •

*The optimal values of metrics are in bold.

Figure 7: Medical claim dataset: the lift plot for GLMs without clustering
method (M2), GLMs with clustering method (M3), mixture model without
clustering method (M5) and mixture model with clustering method (M6) on the
testing data.

##  6 Conclusion

This paper discusses the framework for the implementation of covariate
clustering in mixture regression models. By adding a regularization term for
similar covariates and using the EM-ADMM algorithm to solve the convex
optimization problem, the proposed framework allows us to identify and analyze
the clustering effects of covariates on heterogeneous subgroups of
populations, thus simultaneously exploring the complex network graph of the
covariates and conducting feature selection. The proposed method can also be
considered as a convex clustering method for the simultaneous grouping of the
covariate information and the outcome information. Both simulation and
empirical results show that our model can accurately identify the potential
clusters of covariates and perform better than traditional (or mixture) GLMs
in multiple datasets for the healthcare costs prediction.

Future research can explore the proposed framework from different angles.
Firstly, while the Gamma distribution is utilized as the underlying
distribution in the mixture model to model positive continuous response
variables, it can be extended to other distributional assumptions. For
instance, one may replace the Gamma distribution with Poisson or Tweedie
distributions, to construct predictive models for the number of claims in the
count type or the aggregate claim amounts containing zero values. Secondly,
while the two empirical examples demonstrated the application of the model to
cluster covariates with similar observations and to cluster multiple levels of
a factor-type variable, covariate clustering is also beneficial when reducing
data dimensionality or characterizing the correlation relationship among
covariates. Thus, the proposed model can be extended to more scenarios, such
as spatio-temporal modeling based on various forms of penalties (e.g., well-
known SCAD or MCP penalties).

Future studies can also explore how to enhance the interpretability of the
proposed method, especially when covariate clustering is applied. It would be
of interest to investigate whether the identified clusters can lead to
meaningful insights into the characteristics of the population, or provide
additional information for risk assessment. Furthermore, the proposed model
can be further applied to other fields, such as finance, life insurance, and
environmental science, where the heterogeneity of response variables is also
commonly observed.

## Acknowledgement

Zhengxiao Li acknowledges the financial support from National Natural Science
Fund of China (Grant No. 71901064) and “the Fundamental Research Funds for the
Central Universities” in UIBE (Grant No. CXTD13-02). Yifan Huang acknowledges
the financial support from the National Natural Science Foundation of China
(Grant No. 72201062) and the National Social Science Fundation of China (Grant
No. 22ATJ005).

## Appendices

###  A The objective function in scaled ADMM

At the M-step, we consider that the $w_{h}$ are constant in
$Q(\bm{\Psi};\bm{\Psi}^{(m)})$ and maximize $Q(\bm{\Psi};\bm{\Psi}^{(m)})$
with respect to the other parameters in $\bm{\Psi}$ with the following form:

$\mathrm{maximize}_{\beta_{0,h},\bm{\beta}_{h},\phi_{h}}\quad\sum_{i=1}^{n}\pi_{i,h}^{(m)}\log{f\left(y_{i}|\mu_{h}(\bm{x}_{i}),\phi_{h}\right)}-\omega\gamma{\left\|\bm{\beta}_{h}\right\|_{2}^{2}}-\omega\frac{v}{2}\sum_{j,k\in
I}s_{j,k}\lvert{\beta}_{j,h}-{\beta}_{k,h}\rvert,$ (A.1)

for $h=1,..,H$. Note that in (A.1), each set of parameters
$\\{\beta_{0,h},\bm{\beta}_{h},\phi_{h}\\}$ acts independently on the
objective function, so the optimization problem can be decomposed according to
the components and run in parallel.

For simplicity we omit the subscript $h$ in the following because the updating
process is consistent for $h=1,\dots,H$. The objective can be reformulated as
a maximization problem in the following form:

$\mathrm{maximize}_{\beta_{0},\bm{\beta},\phi}\quad\sum_{i=1}^{n}\pi_{i}^{(m)}\log{f\left(y_{i}|\mu(\bm{x}_{i}),\phi\right)}-\omega\gamma{\left\|\bm{\beta}\right\|_{2}^{2}}-\omega\frac{v}{2}\sum_{j,k\in
I}s_{j,k}\lvert{\beta}_{j}-{\beta}_{k}\rvert,$ (A.2)

summing up to obtain the overall maxima of $Q(\bm{\Psi};\bm{\Psi}^{(m)})$.

Considering the ADMM algorithm is well suited for decomposable dual
optimization problems and has good convergence properties, we introduce a new
set of auxiliary variables $\bm{z}=\\{z_{j,k}\\}$ in (A.2). The aim of M-step
is to take the following form:

$\displaystyle\mathrm{minimize}_{\beta_{0},\bm{\beta},\phi,\bm{z}}\quad$
$\displaystyle-\sum_{i=1}^{n}\pi_{i}^{(m)}\log{f\left(y_{i}|\mu(\bm{x}_{i}),\phi\right)}+\omega\gamma{\left\|\bm{\beta}\right\|_{2}^{2}}+\omega\frac{v}{2}\sum_{j,k\in
I}s_{j,k}\lvert z_{j,k}-z_{k,j}\rvert,$ s.t. $\displaystyle
z_{j,k}-\beta_{j}=0,\quad{\forall}j,k\in\\{1,2,\dots,p\\},$ (A.3)

where $\bm{z}\in\mathbb{R}^{p\times p}$ constitutes a surrogate for the
original regression coefficients
$\bm{\beta}=(\beta_{1},\beta_{2},\dots,\beta_{p})$. The solution to this
problem can be obtained by minimizing the following augmented Lagrangian
function:

$\displaystyle L_{\rho}\left(\beta_{0},\bm{\beta},\phi,\bm{z},\bm{r}\right)=$
$\displaystyle-\sum_{i=1}^{n}\pi_{i}^{(m)}\log{f\left(y_{i}|\mu(\bm{x}_{i}),\phi\right)}+\omega\gamma{\left\|\bm{\beta}\right\|_{2}^{2}}$
$\displaystyle+\omega\frac{v}{2}\sum_{j=1}^{p}\sum_{k=1}^{p}s_{j,k}\lvert
z_{j,k}-z_{k,j}\rvert$
$\displaystyle+\frac{\rho}{2}\sum_{j=1}^{p}\sum_{k=1}^{p}\left(z_{j,k}-\beta_{j}+r_{j,k}\right)^{2}+C,$
(A.4)

where $\bm{r}=\left\\{r_{j,k}\right\\}\in\mathbb{R}^{p\times p}$ is a set of
scaled dual variables, $\rho>0$ is a positive penalty parameter and $C$ is a
constant. Note that (A) is a scaled ADMM, see e.g. Boyd et al. (2011) in
Section 3.1.1. The ADMM algorithm is also widely discussed in recent studies
(Chen et al., 2021, Zhang et al., 2022).

###  B Update the regression coefficients

In (• ‣ 3.1), when considering the commonly used log link function (i.e.,
$\mu_{i}=\mu(\bm{x}_{i})=\exp(\beta_{0}+\sum_{j=1}^{p}x_{i,j}\beta_{j})$), the
derivatives of $g(\cdot)$ based on the chain rule are:

$\displaystyle\frac{\partial g(\beta_{0},\bm{\beta},\phi)}{\partial\beta_{0}}$
$\displaystyle=-\sum_{i=1}^{n}\frac{f_{\mu}^{{}^{\prime}}(y_{i}|\mu_{i},\phi)}{f(y_{i}|\mu_{i},\phi)}\mu_{i},$
(A.5) $\displaystyle\frac{\partial
g(\beta_{0},\bm{\beta},\phi)}{\partial\beta_{j}}$
$\displaystyle=-\sum_{i=1}^{n}\frac{f_{\mu}^{{}^{\prime}}(y_{i}|\mu_{i},\phi)}{f(y_{i}|\mu_{i},\phi)}\cdot\mu_{i}\cdot
x_{i,j}+2\omega\gamma\cdot\beta_{j}$
$\displaystyle-\rho\left[\sum_{k=1}^{p}\left(z_{j,k}^{(m,t)}+r_{j,k}^{(m,t)}\right.)-p\beta_{j}\right],\quad{\forall}j=1,2,\dots,p,$
(A.6)

where $f_{\mu}^{{}^{\prime}}(\cdot)$ represents the derivative function of the
corresponding parametric distribution with respect to its mean parameter. The
Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm can be applied to calculate
the updated regression coefficients $\beta_{0}^{(m,t+1)},\bm{\beta}^{(m,t+1)}$
based on (A.5) and (A.6), and $\phi^{(m,t+1)}$ can be updated by maximum
likelihood estimation in numerical method.

###  C Update the auxiliary variables

According to Hallac et al. (2015), this problem has a closed form solution;
more specifically, for $j,k=1,2,\dots,p$, the updating formulas are:

$\left\\{\begin{aligned}
z_{j,k}^{(m,t+1)}=\theta_{j,k}^{(m,t+1)}\left(\beta_{j}^{(m,t+1)}-r_{j,k}^{(m,t)}\right)+(1-\theta_{j,k}^{(m,t+1)})\left(\beta_{k}^{(m,t+1)}-r_{k,j}^{(m,t)}\right),\\\
z_{k,j}^{(m,t+1)}=(1-\theta_{j,k}^{(m,t+1)})\left(\beta_{j}^{(m,t+1)}-r_{j,k}^{(m,t)}\right)+\theta_{j,k}^{(m,t+1)}\left(\beta_{k}^{(m,t+1)}-r_{k,j}^{(m,t)}\right),\end{aligned}\right.$
(A.7)

with

$\theta_{j,k}^{(m,t+1)}=\begin{cases}&0.5,\quad\text{if}\quad
d_{j,k}^{(m,t+1)}\triangleq\sqrt{\left[(\beta_{j}^{(m,t+1)}-r_{j,k}^{(m,t)})-(\beta_{k}^{(m,t+1)}-r_{k,j}^{(m,t)})\right]^{2}}=0,\\\
&\max{(1-\frac{\omega vs_{j,k}}{\rho
d_{j,k}^{(m,t+1)}},0.5)},\quad\text{else}.\end{cases}$ (A.8)

###  D The description of the questionnaire variables in diabetes expense
dataset

The following part is concerning some question statements for the
questionnaire data, which are related to the knowledge of diabetes. In the V4
statement, it is posited that emotional stress can have an impact on blood
sugar levels in individuals with diabetes. The V5 statement asserts that blood
sugar levels may not necessarily increase the likelihood of developing
complications associated with diabetes. The V77 statement suggests that
smoking can significantly increase the risk of severe diabetic foot amputation
in individuals with diabetes. Finally, the V79 statement notes that the
combination of smoking and diabetes can lead to an increased risk of stroke.
These statements highlight the complex interplay of factors that might affect
the management and potential complications of diabetes, underscoring the
importance of a holistic approach to treatment and risk reduction in
individuals with this condition.

## References

  * Andrade et al. (2021) Daniel Andrade, Kenji Fukumizu, and Yuzuru Okajima.  Convex covariate clustering for classification.  _Pattern Recognition Letters_ , 151:193–199, 2021.
  * Bertsimas et al. (2008) Dimitris Bertsimas, Margrét V Bjarnadóttir, Michael A Kane, J Christian Kryder, Rudra Pandey, Santosh Vempala, and Grant Wang.  Algorithmic prediction of health-care costs.  _Operations Research_ , 56(6):1382–1392, 2008\.
  * Boyd et al. (2011) Stephen Boyd, Neal Parikh, Eric Chu, Borja Peleato, Jonathan Eckstein, et al.  Distributed optimization and statistical learning via the alternating direction method of multipliers.  _Foundations and Trends® in Machine learning_ , 3(1):1–122, 2011.
  * Cao et al. (2012) Qing Cao, Bradley T Ewing, and Mark A Thompson.  Forecasting medical cost inflation rates: A model comparison approach.  _Decision Support Systems_ , 53(1):154–160, 2012\.
  * Chen et al. (2021) Jingxiang Chen, Quoc Tran-Dinh, Michael R Kosorok, and Yufeng Liu.  Identifying heterogeneous effect using latent supervised clustering with adaptive fusion.  _Journal of Computational and Graphical Statistics_ , 30(1):43–54, 2021.
  * Cheng et al. (2022) Chao Cheng, Xingdong Feng, Xiaoguang Li, and Mengyun Wu.  Robust analysis of cancer heterogeneity for high-dimensional data.  _Statistics in Medicine_ , 41(27):5448–5462, 2022\.
  * Chi and Lange (2015) Eric C Chi and Kenneth Lange.  Splitting methods for convex clustering.  _Journal of Computational and Graphical Statistics_ , 24(4):994–1013, 2015.
  * Dempster et al. (1977) Arthur P Dempster, Nan M Laird, and Donald B Rubin.  Maximum likelihood from incomplete data via the em algorithm.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 39(1):1–22, 1977.
  * Devriendt et al. (2021) Sander Devriendt, Katrien Antonio, Tom Reynkens, and Roel Verbelen.  Sparse regression with multi-type regularized feature modeling.  _Insurance: Mathematics and Economics_ , 96:248–261, 2021\.
  * Duncan et al. (2016) Ian Duncan, Michael Loginov, and Michael Ludkovski.  Testing alternative regression frameworks for predictive modeling of health care costs.  _North American Actuarial Journal_ , 20(1):65–87, 2016.
  * Fellingham et al. (2015) Gilbert W Fellingham, Athanasios Kottas, and Brian M Hartman.  Bayesian nonparametric predictive modeling of group health claims.  _Insurance: Mathematics and Economics_ , 60:1–10, 2015\.
  * Girvan and Newman (2002) Michelle Girvan and Mark EJ Newman.  Community structure in social and biological networks.  _Proceedings of the national academy of sciences_ , 99(12):7821–7826, 2002.
  * Gneiting and Raftery (2007) Tilmann Gneiting and Adrian E Raftery.  Strictly proper scoring rules, prediction, and estimation.  _Journal of the American statistical Association_ , 102(477):359–378, 2007.
  * Hallac et al. (2015) David Hallac, Jure Leskovec, and Stephen Boyd.  Network lasso: Clustering and optimization in large graphs.  In _Proceedings of the 21th ACM SIGKDD international conference on knowledge discovery and data mining_ , pages 387–396, 2015.
  * Ickowicz and Sparks (2017) Adrien Ickowicz and Ross Sparks.  Modelling hospital length of stay using convolutive mixtures distributions.  _Statistics in medicine_ , 36(1):122–135, 2017\.
  * Khalili and Chen (2007) Abbas Khalili and Jiahua Chen.  Variable selection in finite mixture of regression models.  _Journal of the american Statistical association_ , 102(479):1025–1038, 2007.
  * Khalili and Lin (2013) Abbas Khalili and Shili Lin.  Regularization in finite mixture of regression models with diverging number of parameters.  _Biometrics_ , 69(2):436–446, 2013.
  * Kurz and Hatfield (2019) Christoph F Kurz and Laura A Hatfield.  Identifying and interpreting subgroups in health care utilization data with count mixture regression models.  _Statistics in medicine_ , 38(22):4423–4435, 2019\.
  * Lee (2021) Simon CK Lee.  Addressing imbalanced insurance data through zero-inflated poisson regression with boosting.  _ASTIN Bulletin: The Journal of the IAA_ , 51(1):27–55, 2021.
  * Ma and Huang (2017) Shujie Ma and Jian Huang.  A concave pairwise fusion approach to subgroup analysis.  _Journal of the American Statistical Association_ , 112(517):410–423, 2017.
  * Mihaylova et al. (2011) Borislava Mihaylova, Andrew Briggs, Anthony O’Hagan, and Simon G Thompson.  Review of statistical methods for analysing healthcare resources and costs.  _Health economics_ , 20(8):897–916, 2011.
  * Miljkovic and Grün (2016) Tatjana Miljkovic and Bettina Grün.  Modeling loss data using mixtures of distributions.  _Insurance: Mathematics and Economics_ , 70:387–396, 2016\.
  * Neelon et al. (2015) Brian Neelon, Fan Li, Lane F Burgette, and Sara E Benjamin Neelon.  A spatiotemporal quantile regression model for emergency department expenditures.  _Statistics in medicine_ , 34(17):2559–2575, 2015\.
  * Shen and He (2015) Juan Shen and Xuming He.  Inference for subgroup analysis with a structured logistic-normal mixture model.  _Journal of the American Statistical Association_ , 110(509):303–312, 2015.
  * Sui et al. (2018) Xiaopeng Lucia Sui, Li Xu, Xiaoning Qian, and Tie Liu.  Convex clustering with metric learning.  _Pattern Recognition_ , 81:575–584, 2018.
  * Tsai and Cheng (2021) Cary Chi-Liang Tsai and Echo Sihan Cheng.  Incorporating statistical clustering methods into mortality models to improve forecasting performances.  _Insurance: Mathematics and Economics_ , 99:42–62, 2021\.
  * Witten et al. (2014) Daniela M Witten, Ali Shojaie, and Fan Zhang.  The cluster elastic net for high-dimensional regression with unknown variable grouping.  _Technometrics_ , 56(1):112–122, 2014.
  * Xie et al. (2016) Yang Xie, Günter Schreier, Michael Hoy, Ying Liu, Sandra Neubauer, David CW Chang, Stephen J Redmond, and Nigel H Lovell.  Analyzing health insurance claims on different timescales to predict days in hospital.  _Journal of biomedical informatics_ , 60:187–196, 2016\.
  * Yin and Lin (2016) Cuihong Yin and X Sheldon Lin.  Efficient estimation of erlang mixtures using iscad penalty with insurance application.  _ASTIN Bulletin: The Journal of the IAA_ , 46(3):779–799, 2016.
  * Zhang et al. (2022) Xin Zhang, Jia Liu, and Zhengyuan Zhu.  Learning coefficient heterogeneity over networks: A distributed spanning-tree-based fused-lasso regression.  _Journal of the American Statistical Association_ , pages 1–13, 2022\.
  * Zhu (2017) Yunzhang Zhu.  An augmented admm algorithm with application to the generalized lasso problem.  _Journal of Computational and Graphical Statistics_ , 26(1):195–204, 2017.