# Anomaly Detection in Health Insurance Claims Using Bayesian Quantile Regression

Authors: Ezekiel N. N. Nortey; Reuben Pometsey; Louis Asiedu; Samuel Iddi; Felix O. Mettle
Journal: International Journal of Mathematics and Mathematical Sciences
 (2021)

---

## Abstract

Research has shown that current health expenditure in most countries, especially in sub-Saharan Africa, is inadequate and unsustainable. Yet, fraud, abuse, and waste in health insurance claims by service providers and subscribers threaten the delivery of quality healthcare. It is therefore imperative to analyze health insurance claim data to identify potentially suspicious claims. Typically, anomaly detection can be posited as a classification problem that requires the use of statistical methods such as mixture models and machine learning approaches to classify data points as either normal or anomalous. Additionally, health insurance claim data are mostly associated with problems of sparsity, heteroscedasticity, multicollinearity, and the presence of missing values. The analyses of such data are best addressed by adopting more robust statistical techniques. In this paper, we utilized the Bayesian quantile regression model to establish the relations between claim outcome of interest and subject-level features and further classify claims as either normal or anomalous. An estimated model component is assumed to inherently capture the behaviors of the response variable. A Bayesian mixture model, assuming a normal mixture of two components, is used to label claims as either normal or anomalous. The model was applied to health insurance data captured on 115 people suffering from various cardiovascular diseases across different states in the USA. Results show that 25 out of 115 claims (21.7%) were potentially suspicious. The overall accuracy of the fitted model was assessed to be 92%. Through the methodological approach and empirical application, we demonstrated that the Bayesian quantile regression is a viable model for anomaly detection.


---

## Body

## 1. Introduction

Providing quality and accessible healthcare is a major political decision that governments make all over the world. A major challenge in achieving this is funding. According to Novignon et al. [1], majority of countries in developing regions, especially sub-Saharan Africa (SSA), rely on donor grants and loans to finance healthcare. Such expenditures are not only unsustainable but also inadequate considering the enormous healthcare burden in the region. Apart from direct government investment into health services, health insurance schemes have been instituted to reduce the proportion of health service cost borne by the populace in many countries. However, fraud, abuse, and waste in health insurance claims by service providers and subscribers threaten the delivery of quality healthcare. The National Healthcare Antifraud Association (NHCAA) of the United States defines healthcare fraud as an intentional deception or misrepresentation made by a person or an entity, with the knowledge that the deception could result in some kinds of unauthorized benefits to that person or entity [2, 3]. This “intentional deception or misrepresentation” can produce observations in the data that do not conform to normal observed patterns. These can therefore be regarded as anomalies or outliers. There are several statistical methods that are employed in the identification of fraud, especially in health insurance claims. Wilson [4] applied the logistic regression model to detect auto insurance fraud. This method is only possible if there is information on both legitimate and fraudulent claims. Moreover, the labeling process may be fraught with error such that any model that is built on this may not automatically capture these errors. Ekina et al. [2] used Bayesian co-clustering model to detect health insurance fraud. This method has been used for the analysis of dyadic data, and since “conspiracy fraud” involves a connection between two entities, the Bayesian co-clustering is an appropriate method. Liu and Vasarhelyi [5] used a clustering method incorporating geo-locating information to detect fraud in Medicare/Medicaid data in the US.Ekin et al. [6] provided a comprehensive assessment of medical fraud research. They mentioned that health fraud research is an emerging field which has gained some attention from researchers but acknowledged that more research needs to be conducted especially using Bayesian techniques. Bayesian methods have generally been proven to be an effective approach of identifying fraud. This is because Bayesian methods do not make strong statistical assumption as the frequentist counterparts.Below are some reasons why Bayesian techniques would be appropriate for health fraud study [6].(i)
Medical data are sparse because fraudulent cases are rare, as compared to legitimate cases.(ii)
Data are also dynamic. Both fraudulent and legitimate patterns evolve over time.(iii)
Fraudulent claims are not homogeneous.(iv)
Missing values are present, and these can lead to over- or undersampling, nonrepresentativeness, and potential bias in inference.(v)
Medical data are skewed and non-normal.To address these challenges, this study proposes the use of Bayesian quantile regression techniques and Bayesian mixture model to identify fraud in health insurance data.This study seeks to specify the posterior distribution of the parameters of interest given the data; identify the distribution of anomalous claims; and estimate the probability of each claim to identify potential anomalous claims.


## 2. Materials and Methods

### 2.1. Bayesian Method

From the Bayesian perspective, parameter vector of a given model is assumed to have a probability distribution referred to as the prior probability distribution. This prior probability distribution is informed by the prior knowledge about the parameter. The data from which the parameter is to be estimated are generated by a probability process, and this generating process is referred to as the likelihood. Of interest is the posterior probability which is the probability of the unknown parameter given the data. GivenYi as a vector of observations, independently and identically distributed, and θi a vector of unknown parameters, then the posterior probability distribution of the parameters is represented as(1)Pθi|Yi=Pθi×PYi|θi∫θiPθi×PYi|θidθi,where Pθi is the prior probability of the parameter θi and Pθi|Yi is the posterior probability of θi given the observed data Yi.Generally, the prior may be chosen based on some information available and hence called informative prior. It may be chosen based on little or no information available. Such priors are called noninformative prior. The uniform distribution or Jeffrey’s prior is mostly used when there is no prior information about the unknown parameter. A prior may also either be a proper prior or an improper prior. A chosen proper prior results in a proper posterior probability. An improper prior is a distribution that does not integrate to one. However, if an improper prior is chosen properly, the resulting posterior distribution may be a proper posterior distribution.There are other ways of choosing the prior. The prior can be chosen such that the resultant posterior distribution has the same form as the prior probability. Such priors are called conjugate.


### 2.2. Quantile Regression

This study used the quantile regression to capture the relationship between the variables and to estimate parameters of the model. This is because of its ability to capture the relationship across the entire data without imposing strong parametric assumptions about the response variable.It is therefore an adequate technique which captures all dynamics in the data and addresses challenges associated with health insurance claims data.Furthermore, given a set of response variablesY and explanatory variables X, the conditional cumulative distribution function is(2)FY|XY|X=x=PY≤y|X=x,and for any τ, 0<τ<1, the τth conditional quantile function is(3)QτY|X=x=infy:FY|Xy≥τ.The quantile functionQτ is monotonic increasing such that for any two quantiles τ1 and τ2 where τ1<τ2, Qτ1<Qτ2.


### 2.3. Description of the Method

To specify the density distribution that describes the anomalous claims, a posterior distribution of the parameters of interest (β, σ, and mi′s) is first constructed using the quantile regression model:(4)Y=Xβ+θmi+wσmiui,where EY=QτY=Xβ+θmi and ui has a standard normal distribution. The parameters to be estimated are the regression coefficients β, a scale parameter σ, and the mixture component mi.From the Bayesian perspective, given thatYi are independently and identically random variables from which the parameters β, σ, and mi′s are to be estimated, the posterior distribution fβ,σ,mi|y1,y2,…,yn is constructed to be proportional to some function gβ,σ,mi expressed as(5)fβ,σ,mi|y1,y2,…,yn∝gβ,σ,mi.Now(6)gβ,σ,mi=fY|β,σ,mi×fβ|Y,σ,mi×fσ|Y,β,mi×fmi|Y,β,σ,by the factorization theorem. To estimate say β, σ and mi are treated as constants and their functions are included as a normalizing constant.Thus,(7)gβ,σ,mi∝fβ|Y,σ,mi,so that(8)fβ,σ,mi|y1,y2,…,yn∝fβ|Y,σ,mi.We then use the Gibbs sampler to sampleβ from the conditional distribution fβ|Y,σ,mi. This is the same as sampling β from the posterior distribution fβ,σ,mi|y1,y2,…,yn. Samples are also drawn from the conditional distributions fσ|Y,β,mi and fmi|Y,β,σ to estimate σ and mi, respectively. It is assumed that the estimates for mi′s form a mixture distribution of unknown densities.


### 2.4. The Model

SupposeYi, i=1,2,…,n are the set of independent variables identically and independently distributed, Xi are the regressors or covariates, and βj are the regression coefficients.The model for the Bayesian quantile regression as proposed by Kozumi and Kobayashi [7] is(9)Yi=Xijβj+εi,i=1,2,…,n,Yi=Xijβj+σθzi+wziui,i=1,2,…,n,where σ is a scale parameter, εi=θσzi+wσziui are the error terms, θ=1−2τ/1−τ, and w2=2/τ1−τ. zi has a standard exponential distribution, and ui has the standard normal distribution. We replace σzi with mi and equation (9) becomes(10)Yi=Xijβj+θmi+wσmiui,which is the same as equation (4).Given a random error variableε with scale parameter σ, asymmetry parameter τ, and location parameter μ, ε is said to have the asymmetry Laplace distribution (assumed by Yu and Moyeed [8]) given by the density(11)fε=τ1−τσexp−ρτεσ,where ρτu=τIu≥0+1−τIu<0u is the loss or check function, τ is on the interval [0, 1], and the τth quantile of ε is zero. The loss function assigns τ of the errors to positive values and 1−τ of the errors to negative values. Now if we let ε=Yi−Xijβj−θmi, then equation (11) becomes(12)fY,X,β,vi=τ1−τσexp−ρτYi−Xijβj−θmiσ.From equation (12), μ=Xijβj+θmi is the location parameter. vi=θmi can only be zero through θ when τ=1/2. The term mi, however, is an exponential random variable and cannot be zero. The presence of anomaly in Y can be demonstrated through mi by assuming that the predicted values of mi have a mixture distribution.It is somewhat intractable to derive the full conditional posterior distributions for the parameters from equation (12). Markov Chain Monte Carlo (MCMC) algorithms such as Gibbs sampler [9] and Metropolis–Hastings algorithm [10] enable one to estimate the parameters from the posterior distribution.MCMC algorithms generate random numbers that will have stationary transition probabilities, hence an equilibrium distribution. In the Bayesian perspective, this equilibrium distribution, given some initial conditions, mimics the posterior distribution. Having done this, we can then compute statistics such as mean and variance from this distribution as estimates of the posterior distribution.Assumingεi has an exponential normal mixture distribution, it can be shown that Yi will be normally distributed with mean, Xijβj+θmi, variance, w2σmi, and a normal likelihood density of the form:(13)LYi=∏i=1nfY|X,β,mi,σ=∏i=1n1wσmiexp−Yi−Xijβj−θmi22w2σmi=exp−∑i=1nYi−Xijβj−θmi22w2σmi∏i=1n1wσmi.The variablesβ,mi, and σ are estimated by assuming the following priors:(14)βj∼NBo,Vo,mi∼exp1σ,σ∼IGno2,so2.By specifying the above conjugate priors for the parameters, their full marginal posterior distributions can be derived and we can subsequently sample from them.From the likelihood function in equation (13) and the above prior distributions, the posterior distribution is given as(15)fβ,mi,σ|Yi∝NXβ+θmi,w2σmi×NBo,Vo×exp1σ×IGno2,so2,which becomes(16)fβ,mi,σ|Yi∝exp−∑i=1nYi−Xijβj−θmi22w2σmi∏i=1n1wσmi×12πVoexp−βj−Bo22Vo×1σnexp−∑i=1nmiσ×1σno/2−1exp−so2σ.The marginal conditional posterior distributions of the parameters are derived from equation (16) as follows:The full conditional posterior distribution ofβ is a normal distribution given as(17)fβ|σ,mi,Yi,Xi∝∏i=1n1σmi×12πVoexp−∑i=1nYi−Xijβj−θmi22w2σmiexp−βj−Bo22Vo=σ−n/2min/22πVoexp−∑i=1nYi−Xijβj−θmi22w2σmi−βj−Bo22Vo=σ−n/2min/22πVoexp−12∑X′Xβ2w2σVo+β2Vo−2XYβw2σmi+2Xβθmiw2σmi−2ββoVo+Kfβ|σ,mi,Yi,Xi∝σ−n/2min/22πVoexp−12∑X′Xw2σmi+1Voβ2−∑XY−θmiw2σmi+βoVoβ+K,where K is a constant that does not depend on β.Note that a functionx↦expax2+bx+c with domain on the whole real line is an unnormalized probability density function (pdf) if and only if a<0. In this case, the unnormalized pdf of the normal distribution with mean (μ) and variance (σ2) is(18)μ=−b2a,σ2=−12a.Comparing equation (17) to x↦expax2+bx+c yields(19)σ2=1∑X′X/w2σmi+1/Vo=∑X′Xw2σmi+1Vo−1,μ=bσ2=∑XY−θmiw2σmi+βoVo∑X′Xw2σmi+1Vo−1.If we letσ2=V¯ and μ=B¯, then the distribution of β is normal with parameters B¯ and V¯ where(20)V¯−1=∑i=1nXi′Xiw2σmi+1Vo,B¯=V¯∑i=1nXiYi−θmiw2σmi+β0Vo.It can also be derived from equation (16) that the full conditional distribution of mi is(21)fmi|β,σ,λi,γ,Yi∝1σn∏1wσmiexp−∑i=1nYi−Xijβj−θmi22w2σmiexp−∑i=1nmiσ∝1min/2exp−∑i=1nY−Xβ2−2Y−Xβθmi−θmi22w2σmi+miσ∝mi−n/2exp−∑i=1n−2Y−Xβθmi2w2θmi+Y−Xβ22w2σmi−θ2mi22w2σmi−miσ∝mi−n/2Kexp∑i=1nY−Xβ22w2σmi−1−∑i=1nθ2w2σ+1σmifmi|β,σ,λi,γ,Yi∝mi−n/2exp−12∑i=1nY−Xβ2w2σmi−1+∑i=1nθw2σ+2σmi.The full conditional posterior distribution of themi is the generalized inverse Gaussian distribution with parameters 1/2, μ^, and ν^ given as(22)fmi|β,σ,Yi,Xi∝mi−n/2exp−12μi^2mi−1+ν^2mi,where μi^2=∑Yi−Xijβj2/σw2 and ν^2=∑2/σ+θ/σw2.Again from equation (16), the full conditional distribution of σ follows the inverse gamma distribution given as(23)fσ|β,mi,Y,X∝1σn/21σn1σno/2−1exp−∑i=1nY−Xβ−θmi22w2σmi−∑i=1nmiσ−so2σ,fσ|β,mi,Y,X∝1σ3n+no/2−1exp−1σso2+∑mi+∑Yi−Xijβj−θmi22w2mi.This is an inverse gamma distribution with shape parameter3n+no/2 and scale parameter so/2+∑mi+∑Yi−Xijβj−θmi2/2w2mi.


### 2.5. Estimating the Mixture Components ofmi

We used theR package for the Bayesian mixture model, “BayesMix,” to estimate the mixture components of mi. BayesMix provides Gibbs sampling of the posterior distribution, a method to set up the model and specify the priors and initial values required for the Gibbs sampler.To run the Gibbs sampler, we specify the prior distributions of the parameters and initial values for the hyperparameters in the proposed models. The regression coefficients are normally distributed with meanβo and variance Vo. We set initial values for these parameters as 0 and 10 (with precision as 1/10), respectively.mi′s are exponentially distributed with rate parameter 1/σ. The prior distribution of the scale parameter σ is an inverse gamma distribution with hyperparameters no/2 and so/2. The quantity no is the prior effective sample size for σ, and so is prior point estimate for σ. For our analysis, we set 1/σ to prec which is the rate of the exponential variable mi. The initial values for prec are generated randomly as uniform values with interval 0 and 1 and 0 and 2. This choice is subjective and arbitrary but with the hope that the data will be informative enough to produce a stationary posterior distribution.


### 2.6. Model Evaluation

The proposed model is evaluated by replicating new values ofY, Yrep and determining how close these values are to the original ones. The approach is called posterior consistency and involves using a posterior predictive distribution.The posterior predictive distribution is given by(24)pYrep|Y=∫ΘfYrep|θi∗fθi|Ydθi,where Θ is our parameter space and θi is the set of parameters we estimated from our model. We sample from this distribution as follows:(1)
Sampleθi from posterior distribution fθi|Y.(2)
Substitute the values ofθi above in fY|θi and draw Yrep.To compare the predictive posterior distribution to our model, we can compute some statistics such as the mean, standard deviation, and order statistics and compare the values. Graphical approaches may also be used. The histogram ofYrep gives the distribution of the predictive posterior distribution, and this can be compared to the model distribution of Y.


#### 2.6.1. Computing the Probability Scores

The densities of the component mixtures are estimated using the Bayesian method. Now, letλ=0 and λ=1 correspond to the density distribution with normal claims and density distribution with anomalous claims (which will have the higher mean), respectively.From the above preamble, we compute the probability of each of the claim belonging to anomalous claim class,PY|λ=1. These probabilities are ranked for all observations, and those with values greater than 0.5 are classified into the anomalous claim class.We treatmi as a mixture of two normal densities and estimate the proportion and density of each component. For μ1 and μ2 (anomalous claims) while the smaller one will be associated with λ=1(anomalous claims) while the smaller one will be associated with λ=0 (normal claims). pY|λi=1, which describes the outlier distribution, can be computed as(25)pY|λ=1=pλ=1pλ=1,Y,σ,mipλ=1pλ=1,Y,σ,mipY,σ,mi+pλ=0pλ=0,Y,σ,mi.


#### 2.6.2. Sensitivity, Specificity, Precision, and Accuracy

The other metrics of performance evaluation can be calculated from Table1.(26)Sensitivity=TPTP+FN=number of true positive classificationsnumber of all positive classifications,specificity=TNFP+TN=number of true negative classificationsnumber of all negative classifications,accuracy=TP+TNTP+FN+FP+TN=number of correct classificationsnumber of all classifications,precision=TPTP+FP.Table 1
Confusion matrix.
ClassificationActualPositiveNegativeTotalPositiveTPFPTP+FPNegativeFNTNFP+TNTotalTP+FNFP+TNTP+TN+FP+FN


### 2.7. Data Acquisition

Secondary data on health insurance claims by 115 people suffering from various cardiovascular diseases located in different states in the USA were used to run the proposed model. The response variable is the claim amount while the covariates/independent variables are location (state), diagnosis, gender, type of care, pharmacy, and type of health facility. All the covariates are categorical variables. Table2 shows the claims data but with the nominal values for the independent variables.Table 2
Claims data with other covariate information.
Claim amountLocationDiagnosisGenderType of carePharmacyType of health facility725711211158352111217023211132928421111281482111329024311143⋮⋮⋮⋮⋮⋮⋮82563111318941211132838531114211731.2321254Location: 1, California; 2, Georgia; 3, Alabama; 4, Virginia. Diagnosis: 1, acute myocardial infarction; 2, heart failure. Gender: 1, female; 2, male. Type of care: 1, outpatient; 2, inpatient. Pharmacy: 1, heparin; 2, nitroglycerin; 3, eptifibatide; 4, clopidogrel; 5, furosemide; 6, metoprolol; 7, bisoprolol; 8, hydrochlorothiazide. Type of health facility: 1, ambulatory service; 2, outpatient clinic; 3, acute care; 4, hospital (source:https://data.cms.gov).


## 2.1. Bayesian Method

From the Bayesian perspective, parameter vector of a given model is assumed to have a probability distribution referred to as the prior probability distribution. This prior probability distribution is informed by the prior knowledge about the parameter. The data from which the parameter is to be estimated are generated by a probability process, and this generating process is referred to as the likelihood. Of interest is the posterior probability which is the probability of the unknown parameter given the data. GivenYi as a vector of observations, independently and identically distributed, and θi a vector of unknown parameters, then the posterior probability distribution of the parameters is represented as(1)Pθi|Yi=Pθi×PYi|θi∫θiPθi×PYi|θidθi,where Pθi is the prior probability of the parameter θi and Pθi|Yi is the posterior probability of θi given the observed data Yi.Generally, the prior may be chosen based on some information available and hence called informative prior. It may be chosen based on little or no information available. Such priors are called noninformative prior. The uniform distribution or Jeffrey’s prior is mostly used when there is no prior information about the unknown parameter. A prior may also either be a proper prior or an improper prior. A chosen proper prior results in a proper posterior probability. An improper prior is a distribution that does not integrate to one. However, if an improper prior is chosen properly, the resulting posterior distribution may be a proper posterior distribution.There are other ways of choosing the prior. The prior can be chosen such that the resultant posterior distribution has the same form as the prior probability. Such priors are called conjugate.


## 2.2. Quantile Regression

This study used the quantile regression to capture the relationship between the variables and to estimate parameters of the model. This is because of its ability to capture the relationship across the entire data without imposing strong parametric assumptions about the response variable.It is therefore an adequate technique which captures all dynamics in the data and addresses challenges associated with health insurance claims data.Furthermore, given a set of response variablesY and explanatory variables X, the conditional cumulative distribution function is(2)FY|XY|X=x=PY≤y|X=x,and for any τ, 0<τ<1, the τth conditional quantile function is(3)QτY|X=x=infy:FY|Xy≥τ.The quantile functionQτ is monotonic increasing such that for any two quantiles τ1 and τ2 where τ1<τ2, Qτ1<Qτ2.


## 2.3. Description of the Method

To specify the density distribution that describes the anomalous claims, a posterior distribution of the parameters of interest (β, σ, and mi′s) is first constructed using the quantile regression model:(4)Y=Xβ+θmi+wσmiui,where EY=QτY=Xβ+θmi and ui has a standard normal distribution. The parameters to be estimated are the regression coefficients β, a scale parameter σ, and the mixture component mi.From the Bayesian perspective, given thatYi are independently and identically random variables from which the parameters β, σ, and mi′s are to be estimated, the posterior distribution fβ,σ,mi|y1,y2,…,yn is constructed to be proportional to some function gβ,σ,mi expressed as(5)fβ,σ,mi|y1,y2,…,yn∝gβ,σ,mi.Now(6)gβ,σ,mi=fY|β,σ,mi×fβ|Y,σ,mi×fσ|Y,β,mi×fmi|Y,β,σ,by the factorization theorem. To estimate say β, σ and mi are treated as constants and their functions are included as a normalizing constant.Thus,(7)gβ,σ,mi∝fβ|Y,σ,mi,so that(8)fβ,σ,mi|y1,y2,…,yn∝fβ|Y,σ,mi.We then use the Gibbs sampler to sampleβ from the conditional distribution fβ|Y,σ,mi. This is the same as sampling β from the posterior distribution fβ,σ,mi|y1,y2,…,yn. Samples are also drawn from the conditional distributions fσ|Y,β,mi and fmi|Y,β,σ to estimate σ and mi, respectively. It is assumed that the estimates for mi′s form a mixture distribution of unknown densities.


## 2.4. The Model

SupposeYi, i=1,2,…,n are the set of independent variables identically and independently distributed, Xi are the regressors or covariates, and βj are the regression coefficients.The model for the Bayesian quantile regression as proposed by Kozumi and Kobayashi [7] is(9)Yi=Xijβj+εi,i=1,2,…,n,Yi=Xijβj+σθzi+wziui,i=1,2,…,n,where σ is a scale parameter, εi=θσzi+wσziui are the error terms, θ=1−2τ/1−τ, and w2=2/τ1−τ. zi has a standard exponential distribution, and ui has the standard normal distribution. We replace σzi with mi and equation (9) becomes(10)Yi=Xijβj+θmi+wσmiui,which is the same as equation (4).Given a random error variableε with scale parameter σ, asymmetry parameter τ, and location parameter μ, ε is said to have the asymmetry Laplace distribution (assumed by Yu and Moyeed [8]) given by the density(11)fε=τ1−τσexp−ρτεσ,where ρτu=τIu≥0+1−τIu<0u is the loss or check function, τ is on the interval [0, 1], and the τth quantile of ε is zero. The loss function assigns τ of the errors to positive values and 1−τ of the errors to negative values. Now if we let ε=Yi−Xijβj−θmi, then equation (11) becomes(12)fY,X,β,vi=τ1−τσexp−ρτYi−Xijβj−θmiσ.From equation (12), μ=Xijβj+θmi is the location parameter. vi=θmi can only be zero through θ when τ=1/2. The term mi, however, is an exponential random variable and cannot be zero. The presence of anomaly in Y can be demonstrated through mi by assuming that the predicted values of mi have a mixture distribution.It is somewhat intractable to derive the full conditional posterior distributions for the parameters from equation (12). Markov Chain Monte Carlo (MCMC) algorithms such as Gibbs sampler [9] and Metropolis–Hastings algorithm [10] enable one to estimate the parameters from the posterior distribution.MCMC algorithms generate random numbers that will have stationary transition probabilities, hence an equilibrium distribution. In the Bayesian perspective, this equilibrium distribution, given some initial conditions, mimics the posterior distribution. Having done this, we can then compute statistics such as mean and variance from this distribution as estimates of the posterior distribution.Assumingεi has an exponential normal mixture distribution, it can be shown that Yi will be normally distributed with mean, Xijβj+θmi, variance, w2σmi, and a normal likelihood density of the form:(13)LYi=∏i=1nfY|X,β,mi,σ=∏i=1n1wσmiexp−Yi−Xijβj−θmi22w2σmi=exp−∑i=1nYi−Xijβj−θmi22w2σmi∏i=1n1wσmi.The variablesβ,mi, and σ are estimated by assuming the following priors:(14)βj∼NBo,Vo,mi∼exp1σ,σ∼IGno2,so2.By specifying the above conjugate priors for the parameters, their full marginal posterior distributions can be derived and we can subsequently sample from them.From the likelihood function in equation (13) and the above prior distributions, the posterior distribution is given as(15)fβ,mi,σ|Yi∝NXβ+θmi,w2σmi×NBo,Vo×exp1σ×IGno2,so2,which becomes(16)fβ,mi,σ|Yi∝exp−∑i=1nYi−Xijβj−θmi22w2σmi∏i=1n1wσmi×12πVoexp−βj−Bo22Vo×1σnexp−∑i=1nmiσ×1σno/2−1exp−so2σ.The marginal conditional posterior distributions of the parameters are derived from equation (16) as follows:The full conditional posterior distribution ofβ is a normal distribution given as(17)fβ|σ,mi,Yi,Xi∝∏i=1n1σmi×12πVoexp−∑i=1nYi−Xijβj−θmi22w2σmiexp−βj−Bo22Vo=σ−n/2min/22πVoexp−∑i=1nYi−Xijβj−θmi22w2σmi−βj−Bo22Vo=σ−n/2min/22πVoexp−12∑X′Xβ2w2σVo+β2Vo−2XYβw2σmi+2Xβθmiw2σmi−2ββoVo+Kfβ|σ,mi,Yi,Xi∝σ−n/2min/22πVoexp−12∑X′Xw2σmi+1Voβ2−∑XY−θmiw2σmi+βoVoβ+K,where K is a constant that does not depend on β.Note that a functionx↦expax2+bx+c with domain on the whole real line is an unnormalized probability density function (pdf) if and only if a<0. In this case, the unnormalized pdf of the normal distribution with mean (μ) and variance (σ2) is(18)μ=−b2a,σ2=−12a.Comparing equation (17) to x↦expax2+bx+c yields(19)σ2=1∑X′X/w2σmi+1/Vo=∑X′Xw2σmi+1Vo−1,μ=bσ2=∑XY−θmiw2σmi+βoVo∑X′Xw2σmi+1Vo−1.If we letσ2=V¯ and μ=B¯, then the distribution of β is normal with parameters B¯ and V¯ where(20)V¯−1=∑i=1nXi′Xiw2σmi+1Vo,B¯=V¯∑i=1nXiYi−θmiw2σmi+β0Vo.It can also be derived from equation (16) that the full conditional distribution of mi is(21)fmi|β,σ,λi,γ,Yi∝1σn∏1wσmiexp−∑i=1nYi−Xijβj−θmi22w2σmiexp−∑i=1nmiσ∝1min/2exp−∑i=1nY−Xβ2−2Y−Xβθmi−θmi22w2σmi+miσ∝mi−n/2exp−∑i=1n−2Y−Xβθmi2w2θmi+Y−Xβ22w2σmi−θ2mi22w2σmi−miσ∝mi−n/2Kexp∑i=1nY−Xβ22w2σmi−1−∑i=1nθ2w2σ+1σmifmi|β,σ,λi,γ,Yi∝mi−n/2exp−12∑i=1nY−Xβ2w2σmi−1+∑i=1nθw2σ+2σmi.The full conditional posterior distribution of themi is the generalized inverse Gaussian distribution with parameters 1/2, μ^, and ν^ given as(22)fmi|β,σ,Yi,Xi∝mi−n/2exp−12μi^2mi−1+ν^2mi,where μi^2=∑Yi−Xijβj2/σw2 and ν^2=∑2/σ+θ/σw2.Again from equation (16), the full conditional distribution of σ follows the inverse gamma distribution given as(23)fσ|β,mi,Y,X∝1σn/21σn1σno/2−1exp−∑i=1nY−Xβ−θmi22w2σmi−∑i=1nmiσ−so2σ,fσ|β,mi,Y,X∝1σ3n+no/2−1exp−1σso2+∑mi+∑Yi−Xijβj−θmi22w2mi.This is an inverse gamma distribution with shape parameter3n+no/2 and scale parameter so/2+∑mi+∑Yi−Xijβj−θmi2/2w2mi.


## 2.5. Estimating the Mixture Components ofmi

We used theR package for the Bayesian mixture model, “BayesMix,” to estimate the mixture components of mi. BayesMix provides Gibbs sampling of the posterior distribution, a method to set up the model and specify the priors and initial values required for the Gibbs sampler.To run the Gibbs sampler, we specify the prior distributions of the parameters and initial values for the hyperparameters in the proposed models. The regression coefficients are normally distributed with meanβo and variance Vo. We set initial values for these parameters as 0 and 10 (with precision as 1/10), respectively.mi′s are exponentially distributed with rate parameter 1/σ. The prior distribution of the scale parameter σ is an inverse gamma distribution with hyperparameters no/2 and so/2. The quantity no is the prior effective sample size for σ, and so is prior point estimate for σ. For our analysis, we set 1/σ to prec which is the rate of the exponential variable mi. The initial values for prec are generated randomly as uniform values with interval 0 and 1 and 0 and 2. This choice is subjective and arbitrary but with the hope that the data will be informative enough to produce a stationary posterior distribution.


## 2.6. Model Evaluation

The proposed model is evaluated by replicating new values ofY, Yrep and determining how close these values are to the original ones. The approach is called posterior consistency and involves using a posterior predictive distribution.The posterior predictive distribution is given by(24)pYrep|Y=∫ΘfYrep|θi∗fθi|Ydθi,where Θ is our parameter space and θi is the set of parameters we estimated from our model. We sample from this distribution as follows:(1)
Sampleθi from posterior distribution fθi|Y.(2)
Substitute the values ofθi above in fY|θi and draw Yrep.To compare the predictive posterior distribution to our model, we can compute some statistics such as the mean, standard deviation, and order statistics and compare the values. Graphical approaches may also be used. The histogram ofYrep gives the distribution of the predictive posterior distribution, and this can be compared to the model distribution of Y.


### 2.6.1. Computing the Probability Scores

The densities of the component mixtures are estimated using the Bayesian method. Now, letλ=0 and λ=1 correspond to the density distribution with normal claims and density distribution with anomalous claims (which will have the higher mean), respectively.From the above preamble, we compute the probability of each of the claim belonging to anomalous claim class,PY|λ=1. These probabilities are ranked for all observations, and those with values greater than 0.5 are classified into the anomalous claim class.We treatmi as a mixture of two normal densities and estimate the proportion and density of each component. For μ1 and μ2 (anomalous claims) while the smaller one will be associated with λ=1(anomalous claims) while the smaller one will be associated with λ=0 (normal claims). pY|λi=1, which describes the outlier distribution, can be computed as(25)pY|λ=1=pλ=1pλ=1,Y,σ,mipλ=1pλ=1,Y,σ,mipY,σ,mi+pλ=0pλ=0,Y,σ,mi.


### 2.6.2. Sensitivity, Specificity, Precision, and Accuracy

The other metrics of performance evaluation can be calculated from Table1.(26)Sensitivity=TPTP+FN=number of true positive classificationsnumber of all positive classifications,specificity=TNFP+TN=number of true negative classificationsnumber of all negative classifications,accuracy=TP+TNTP+FN+FP+TN=number of correct classificationsnumber of all classifications,precision=TPTP+FP.Table 1
Confusion matrix.
ClassificationActualPositiveNegativeTotalPositiveTPFPTP+FPNegativeFNTNFP+TNTotalTP+FNFP+TNTP+TN+FP+FN


## 2.6.1. Computing the Probability Scores

The densities of the component mixtures are estimated using the Bayesian method. Now, letλ=0 and λ=1 correspond to the density distribution with normal claims and density distribution with anomalous claims (which will have the higher mean), respectively.From the above preamble, we compute the probability of each of the claim belonging to anomalous claim class,PY|λ=1. These probabilities are ranked for all observations, and those with values greater than 0.5 are classified into the anomalous claim class.We treatmi as a mixture of two normal densities and estimate the proportion and density of each component. For μ1 and μ2 (anomalous claims) while the smaller one will be associated with λ=1(anomalous claims) while the smaller one will be associated with λ=0 (normal claims). pY|λi=1, which describes the outlier distribution, can be computed as(25)pY|λ=1=pλ=1pλ=1,Y,σ,mipλ=1pλ=1,Y,σ,mipY,σ,mi+pλ=0pλ=0,Y,σ,mi.


## 2.6.2. Sensitivity, Specificity, Precision, and Accuracy

The other metrics of performance evaluation can be calculated from Table1.(26)Sensitivity=TPTP+FN=number of true positive classificationsnumber of all positive classifications,specificity=TNFP+TN=number of true negative classificationsnumber of all negative classifications,accuracy=TP+TNTP+FN+FP+TN=number of correct classificationsnumber of all classifications,precision=TPTP+FP.Table 1
Confusion matrix.
ClassificationActualPositiveNegativeTotalPositiveTPFPTP+FPNegativeFNTNFP+TNTotalTP+FNFP+TNTP+TN+FP+FN


## 2.7. Data Acquisition

Secondary data on health insurance claims by 115 people suffering from various cardiovascular diseases located in different states in the USA were used to run the proposed model. The response variable is the claim amount while the covariates/independent variables are location (state), diagnosis, gender, type of care, pharmacy, and type of health facility. All the covariates are categorical variables. Table2 shows the claims data but with the nominal values for the independent variables.Table 2
Claims data with other covariate information.
Claim amountLocationDiagnosisGenderType of carePharmacyType of health facility725711211158352111217023211132928421111281482111329024311143⋮⋮⋮⋮⋮⋮⋮82563111318941211132838531114211731.2321254Location: 1, California; 2, Georgia; 3, Alabama; 4, Virginia. Diagnosis: 1, acute myocardial infarction; 2, heart failure. Gender: 1, female; 2, male. Type of care: 1, outpatient; 2, inpatient. Pharmacy: 1, heparin; 2, nitroglycerin; 3, eptifibatide; 4, clopidogrel; 5, furosemide; 6, metoprolol; 7, bisoprolol; 8, hydrochlorothiazide. Type of health facility: 1, ambulatory service; 2, outpatient clinic; 3, acute care; 4, hospital (source:https://data.cms.gov).


## 3. Results and Discussion

We run the Gibbs sampler at quantile values,τ=0.95,0.75,0.50,0.25, to observe the trend of the changes in β′s (regression coefficients) and σ (scale parameter). Table 3 shows the estimates of the regression coefficients at the specified quantiles (τ=0.95,0.75,0.50,0.25).Table 3
Posterior parameter estimate at various quantile values.
Parameterτ=0.95τ=0.75τ=0.5τ=0.25β18.8324568.810728.5868956.582426β20.0444650.0562380.0100740.09637β30.4167210.273557−0.035930.23371β40.0946410.127956−0.001630.229348β50.0309740.2178730.3462220.393446β6−0.03138−0.03801−0.0054−0.06249β7−0.001340.000169−0.008190.009342σ8.9801033.0623483.3137970.803441From Table3, the quantile value τ=0.75 is preferred to estimate the anomalous claims in the data. The values of the estimated parameters at τ=0.75 help satisfy the assumption that an anomalous/fraudulent health insurance claim is usually located at the right of a legitimate claim. This suggests that we can estimate the parameters of the proposed model even at τ=0.95. However, the potential reduction factor at τ=0.75 and τ=0.95 is 1.01 and 1.08, respectively. This makes τ=0.75 the preferred quantile value.The mean, standard deviation, and standard error of the mean for each parameter estimated atτ=0.75 are given in Table 4.Table 4
Empirical mean and standard deviation for each parameter, with standard error of the mean atτ=0.75.
ParameterMeanSDNaive SETime-series SEβ[1]8.81072040.249250.0005570.005705β[2]0.05623770.060560.0001350.000989β[3]0.27355690.789520.0017650.084209β[4]0.12795550.103450.0002310.001989β[5]0.21787290.783660.0017520.075442β[6]−0.03801310.049040.000110.000912β[7]0.00016920.084090.0001880.001818mi[1]3.55367951.65510.0037010.008894mi[2]3.83910551.758220.0039320.008311mi[3]3.24036071.481880.0033140.005923mi[4]2.80098531.296630.0028990.005247mi[5]2.90375411.336270.0029880.005209⋮⋮⋮⋮⋮σ3.062346326980.000730.00568We provide the residual and normal Q-Q plots of the model in Figure1 and also the trace plot and density of β and σ in Figure 2.Figure 1
Residuals and normal Q-Q plots.
(a)(b)Figure 2
Trace plot and density ofβ and σ at τ=0.75.The results from Figure2 (diagram of densities) satisfy the assumption of inverse gamma and normal prior distributions made for σ and β, respectively.Our objective is to identify potential fraudulent claims through the estimated values of the unobserved variablemi. Figure 3 shows the histogram and density of mi at τ=0.75.Figure 3
Density ofmi.


### 3.1. Identifying Anomalous Claims

The “BayesMix package” inR was used to set up an MCMC sampler to identify the two components of the mixture of mi′s. The mean (μ), sigma (σ), and proportion (γ) of each component were estimated. Note that the probability of a claim being anomalous is Pλ=1=γ2 and the probability of a claim being normal is Pλ=0=γ1. The results are shown in Table 5.Table 5
Empirical mean, standard deviation, and 95%μ, γ, and σ2.
ParameterMeanSD2.5%97.5%μ110.480.62729.13811.08μ214.172.037910.96418.73γ10.83490.250.0290280.9926γ20.16510.250.0073770.9710σ23.42570.80052.20584.8034The trace plots and densities of the mixture parametersμ and σ2 are shown in Figures 4 and 5, respectively.Figure 4
Trace plot and density ofμ.Figure 5
Trace plot and density ofσ2.The probability of each of the claims belonging to anomalous class,PY|λ=1, is through mi′s, so we use μ2=14.17 for density of PY,m,λ=1 and Pλ=1=0.1651. Moreover, the density of Pλ=0PY,m,λ=0 is μ1=10.48 and Pλ=0=0.8349. Assuming the mixture components are normally distributed, we compute PY|λ=1 using equation (25) and provide ranked results with claims in Table 6. It is evident from Table 6 that each of the first 10 claims with claim amounts 9440,9638,9316,…,9855 has PY|λ=1 greater than 0.5; hence, they would each be classified into the anomalous class.Table 6
Claims and their probabilities of being anomalous.
miClaim amountpY,λ=0,mipY,λ=1,mipY|λ=1mi[36]94402.67E − 150.0002331mi[103]96383.17E − 110.000631mi[45]93161.24E − 100.0007350.999999mi[8]89412.95E − 090.0010630.999982mi[89]91173.30E − 090.0010780.999981mi[83]92413.65E − 090.0010910.999979mi[6]90241.55E − 080.0013010.999924mi[12]99374.02E − 080.0014650.999826mi[4]92848.22E − 080.0016040.999675mi[92]98552.92E − 070.001890.999019⋮⋮⋮⋮⋮mi[41]66610.4072140.0301850.011533mi[93]66240.9355360.0663750.011045mi[82]122720.4669050.0324710.010828mi[55]73560.5232070.0347260.010339mi[75]9562.80.8851950.058080.010222mi[90]61150.8817620.0576460.010186mi[80]79390.8784150.0572340.010152mi[1]72570.5995790.0380330.009886mi[24]65220.626090.0392740.009777mi[70]9712.30.8113450.0505880.009719mi[16]74710.6516350.0405280.009695mi[61]10472.80.803520.0499540.009691mi[50]72060.6753060.041750.009638mi[97]11065.60.7813460.048270.009631mi[47]68310.6878360.0424240.009615The result of our analysis using classification metrics is summarised in Table7.Table 7
Performance metrics.
MetricsValue (%)Sensitivity100Specificity90Accuracy92Precision74


## 3.1. Identifying Anomalous Claims

The “BayesMix package” inR was used to set up an MCMC sampler to identify the two components of the mixture of mi′s. The mean (μ), sigma (σ), and proportion (γ) of each component were estimated. Note that the probability of a claim being anomalous is Pλ=1=γ2 and the probability of a claim being normal is Pλ=0=γ1. The results are shown in Table 5.Table 5
Empirical mean, standard deviation, and 95%μ, γ, and σ2.
ParameterMeanSD2.5%97.5%μ110.480.62729.13811.08μ214.172.037910.96418.73γ10.83490.250.0290280.9926γ20.16510.250.0073770.9710σ23.42570.80052.20584.8034The trace plots and densities of the mixture parametersμ and σ2 are shown in Figures 4 and 5, respectively.Figure 4
Trace plot and density ofμ.Figure 5
Trace plot and density ofσ2.The probability of each of the claims belonging to anomalous class,PY|λ=1, is through mi′s, so we use μ2=14.17 for density of PY,m,λ=1 and Pλ=1=0.1651. Moreover, the density of Pλ=0PY,m,λ=0 is μ1=10.48 and Pλ=0=0.8349. Assuming the mixture components are normally distributed, we compute PY|λ=1 using equation (25) and provide ranked results with claims in Table 6. It is evident from Table 6 that each of the first 10 claims with claim amounts 9440,9638,9316,…,9855 has PY|λ=1 greater than 0.5; hence, they would each be classified into the anomalous class.Table 6
Claims and their probabilities of being anomalous.
miClaim amountpY,λ=0,mipY,λ=1,mipY|λ=1mi[36]94402.67E − 150.0002331mi[103]96383.17E − 110.000631mi[45]93161.24E − 100.0007350.999999mi[8]89412.95E − 090.0010630.999982mi[89]91173.30E − 090.0010780.999981mi[83]92413.65E − 090.0010910.999979mi[6]90241.55E − 080.0013010.999924mi[12]99374.02E − 080.0014650.999826mi[4]92848.22E − 080.0016040.999675mi[92]98552.92E − 070.001890.999019⋮⋮⋮⋮⋮mi[41]66610.4072140.0301850.011533mi[93]66240.9355360.0663750.011045mi[82]122720.4669050.0324710.010828mi[55]73560.5232070.0347260.010339mi[75]9562.80.8851950.058080.010222mi[90]61150.8817620.0576460.010186mi[80]79390.8784150.0572340.010152mi[1]72570.5995790.0380330.009886mi[24]65220.626090.0392740.009777mi[70]9712.30.8113450.0505880.009719mi[16]74710.6516350.0405280.009695mi[61]10472.80.803520.0499540.009691mi[50]72060.6753060.041750.009638mi[97]11065.60.7813460.048270.009631mi[47]68310.6878360.0424240.009615The result of our analysis using classification metrics is summarised in Table7.Table 7
Performance metrics.
MetricsValue (%)Sensitivity100Specificity90Accuracy92Precision74


## 4. Conclusions and Recommendations

The model identified 25 claims (about 22%) estimated at quantile levelτ=0.75 as belonging to the anomalous class and can thus be assumed to be fraudulent claims. The first 10 of these are shown in Table 6.The model correctly identified all anomalous claims (sensitivity of100%) and gave a high specificity value of 90%. The overall accuracy of the model was 92% with a precision of 74%. This study therefore recommends Bayesian quantile regression as a viable model for anomaly detection. It is worthy of note that the study assumed a normal density mixture for the mi′s to be able to estimate the anomalous cluster in the data. Therefore, a choice of a different distribution for the mi′s could result in different findings. Also, instead of using a mixture model, a variable selection method through data augmentation could be used to aid in selecting only appropriate variables for the model.



---

*Source: 6667671-2021-02-24.xml*