Bayesian LASSO with Categorical Predictors: Coding Strategies, Uncertainty Quantification, and Healthcare Applications

1. Introduction

Categorical predictors are widely used in healthcare data analysis as many key variables in medical research are naturally categorical, such as disease status, treatment groups, and demographic factors like gender, race, and socioeconomic status. These predictors allow researchers to assess differences across groups, identify risk factors, and tailor medical treatments to specific populations. Additionally, categorical variables often play a crucial role in clinical decision-making, where classifications like disease severity or diagnostic test results influence patient management. Properly incorporating categorical predictors in statistical models enables robust analysis, guiding evidence-based healthcare policies and personalized medicine. Different from numerical predictors, categorical predictors usually need to be transformed into a numeric format through an appropriate encoding method before being included in the statistical models ([2]). This transformation is necessary because statistical models typically require numerical inputs for computation. It must be decided how to encode categorical variables, with common methods including dummy coding, Helmert coding, and orthogonal contrasts, each emphasizing different aspects of category comparisons. The choice of encoding affects the interpretation of model coefficients, as different schemes highlight various contrasts among categories, such as comparisons to a reference group or overall mean differences. Consequently, a single categorical predictor can be represented in multiple ways within a model, influencing the insights drawn from the analysis and requiring careful selection to align with the study’s objectives.

The impact of coding strategy on variable selection and prediction with linear regression and least absolute shrinkage and selection operator, or LASSO ([3]), has been investigated in [1], which shows that although different coding schemes for categorical predictors do not affect the performance of linear regression, they do impact the performance of LASSO. Despite the success of their study, the results are based on using the point estimates produced from LASSO regression. Uncertainty quantification and statistical inference of LASSO regression using categorical predictors under different coding systems have not been considered. Such a technical gap has motivated us to tackle the problem from an Bayesian perspective, as it is well-acknowledged that fully Bayesian analysis can yield the entire posterior distributions of model parameters (including regression coefficients), therefore statistical inference with Bayesian credible intervals, posterior inclusion probabilities ([4]), and hypothesis testing can be readily performed for a more comprehensive uncertainty assessment ([5]). Specifically, we propose to assess the performance of Bayesian LASSO ([6]) in regression analysis with categorical predictors, which complements the analysis in [1] by incorporating inference procedures. Furthermore, our limited literature search reveals that Bayesian analysis and inference with categorical predictors have rarely been reported ([7]).

Beyond its methodological contributions, our proposed analysis is also novel in terms of the data utilized. In particular, Bayesian techniques have been seldom applied in analyzing chronic autoimmune diseases such as Multiple Sclerosis (MS) ([8,9]). Multiple Sclerosis (MS) is a complex and chronic autoimmune neuroinflammatory disorder that affects the central nervous system and leads to a wide range of physical, cognitive, and emotional impairments. The disease typically manifests in early adulthood and progresses over time, contributing to long-term disability and substantially diminishing health-related quality of life (HRQoL). MS also places a significant economic burden on patients, families, and healthcare systems due to ongoing medical treatments, hospitalizations, and loss of productivity ([10,11]). Understanding the predictors and risk factors associated with MS is essential for early detection, personalized management, and resource allocation. However, many of the relevant predictors—such as sex, race, region, educational attainment, and degree status—are categorical variables, requiring careful statistical handling. Exploring the performance and inference with Bayesian approaches will lead to more robust, interpretable, and clinically meaningful insights into the factors that influence MS onset, progression, and outcomes.

In this study, we investigate the impact of different coding strategies on categorical predictors in prediction, estimation, and inference procedures with Bayesian LASSO using the 2017–2022 Medical Expenditure Panel Survey (MEPS) from the Household Component (HC) and the corresponding Full-Year Consolidated data files. MEPS is a nationally representative survey of the U.S. civilian noninstitutionalized population, collecting detailed information on healthcare utilization, expenditures, insurance coverage, and sociodemographic characteristics. The study sample includes adult individuals aged 18 years and older, with and without a diagnosis of MS. Adults without MS were included in the non-MS group, which allows for the examination of sociodemographic characteristics, health conditions, and healthcare access between MS and non-MS adult populations in the US. Health-related quality of life (HRQoL) was measured using the Veterans RAND 12-Item Health Survey (VR-12) in the MEPS-HC, which yields two standardized scores: the Physical Component Summary (PCS) and the Mental Component Summary (MCS) ([12]). In this study, the response variable of interest is the PCS score, which reflects general physical health, activity limitations, role limitations due to physical health and pain. The dataset includes categorical predictors like sex, race, region, education degree, insurance level, age, and Number of Elixhauser comorbidity conditions ([13]), among which there are categorical predictors such as sex and education, as well as continuous predictors such as age and exhauster Number of Elixhauser comorbidity conditions.

The structure of the article is as follows. First, we demonstrate four classic coding strategies, dummy coding, deviation coding, sequential coding, and Helmert coding, through fitting Bayesian linear regression with simulated data. We then provide a brief introduction to linear regression and LASSO, and use the inference procedure to motivate Bayesian LASSO. In a case study of MEPS data, Bayesian Lasso and alternative methods including linear regression have been to the data for estimation, prediction, and statistical inference. Bayesian analysis yields promising numeric results that have important practical implications, making it particularly powerful in uncertainty quantification, which provides insights for decision-making and scientific interpretation in complex disease analysis.

2. Bayesian Linear Regression with Categorical Predictors

2.1. Coding Strategies with Categorical Predictors

In linear regression analysis, dummy coding is a routine practice to convert the categorical predictor into a group of binary indicators (or dummy variables) which are coded based on choosing a baseline category a priori. In simple linear regression models with standard binary coding, the regression coefficients of binary indicators represent the mean difference in the response variable between the two groups. Therefore, testing the significance of the regression coefficients indicates whether the difference between the two categories is statistically significant. These important statistical implications could potentially benefit healthcare analysis where linear regression with categorical predictors is widely used to include variables such as disease status, treatment groups, and risk classifications. Unlike continuous variables, categorical variables do not have an inherent numerical structure, making direct inclusion infeasible. Encoding strategies, such as dummy coding, transform categorical predictors into binary or numerical representations that capture differences among categories or relationships within the data. Therefore, it improves the interpretability and statistical implications of the corresponding regression coefficients. Additionally, the choice of encoding strategy influences how group comparisons are made, impacting hypothesis testing and inference. Properly encoding categorical predictors is essential for ensuring accurate estimation, meaningful statistical interpretation, and reliable predictions in regression models. Here we adopted four coding strategies for comparison, which are dummy, deviation, sequential and Helmert coding ([14,15]). Regarding the dummy coding we have just discussed, it creates $k-1$ binary indicators for a categorical variable with k categories, where 0’s indicate the baseline category and 1’s indicate one of the remaining $k-1$ categories. Therefore, it is feasible to assess group differences with respect to the baseline level. In addition to the dummy coding that has been widely adopted in regression analysis, other types of coding schemes have also been adopted in areas such as experimental psychology, educational measurement and biomedical studies, where they are used to facilitate specific hypothesis testing, compare ordered categories, or interpret effects relative to the grand mean ([14,16]). For example, deviation coding contrasts each category with the overall mean, providing insight into how each category deviates from the average. It is the same as dummy coding except that the reference group is coded with -1 instead of 0. Besides, sequential coding compares each category to the preceding (or following) one, making it particularly useful for modeling predictors with ordered categories. Specifically, the reference category is coded with 0’s across all binary indicators, and each subsequent category adds a ’1’ to one additional indicator compared to the previous category. Furthermore, Helmert coding compares each category to the mean of all subsequent categories, allowing for progressive contrasts in hierarchical analyses. For all the above coding strategies, a categorical variable with k categories leads to $k-1$ indicators.

To better illustrate the idea, we apply the above 4 coding strategies to encode a categorical variable Education with 5 categories (no degree, high school or GED, bachelor’s degree, master’s and doctorate degree, other degree), from MEPS data where the Physical Component Summary scores, or PCS score, is used as the response variable. These four coding strategies are given in the Table 1, Table 2, Table 3, Table 4, respectively. The five categories of the predictor Education Level are converted into four indicators. The interpretation of the associated coefficients depends on the coding strategy employed. For example, if dummy coding is used by setting “No Degree” as the reference group, the coefficient of binary indicator $D_3$ in Table 1 represents the difference in the average PCS score between individuals with a master’s or doctorate degree and those without a degree. Under deviation coding, the coefficient of $D_3$ in Table 2 indicates how the PCS score of individuals with a bachelor’s degree deviates from the overall mean of all educational levels. In sequential coding, each coefficient represents the contrast between a given category and the category immediately preceding it in the ordered hierarchy. As shown in Table 3, the binary indicator $S_3$ takes the value 1 for individuals with a Master’s or Doctorate degree, as well as those in any subsequent category (e.g., “Other Degree”), and 0 for all lower categories. Consequently, the coefficient for $S_3$ quantifies the mean difference in the outcome between individuals whose highest educational attainment is at least a Master’s or Doctorate degree and those whose highest degree is a Bachelor’s degree. In Helmert coding, each coefficient represents the contrast between a given category and the mean of all subsequent higher categories in the hierarchy. As shown in Table 4, the binary contrast $H_3$ compares individuals with a Bachelor’s degree to the average of those with a Master’s or Doctorate degree and those with another type of degree. Accordingly, the coefficient for $H_3$ quantifies the mean difference in the PCS score between these two groups, indicating how much higher or lower the score is for Bachelor’s degree holders relative to the average score of individuals in the higher education categories. Different coding strategies generate different predictor variables, and the interpretation of each coefficient is strictly determined by the specific coding scheme applied in the regression model.

2.2. Investigation of the Coding Strategies with Bayesian Linear Regression

Reference [1] have demonstrated how to interpret estimated coefficients from linear regression using different coding schemes. Here we provide the Bayesian alternative through the following working example with a continuous response variable y and a categorical predictor of five levels. Unlike the standard linear regression that yields the point estimates of least square regression coefficients, Bayesian linear regression returns the entire posterior distribution of regression coefficients from Markov Chain Monte Carlo (MCMC) so exact statistical inference can be performed even on finite samples.

Using dummy coding, we convert the predictors with five categories group 1 to group 5 into four dummy variables, which can be used to simulate y through the following model:

$$y={\beta }_0+{\beta }_1{x}_1+{\beta }_2{x}_2+{\beta }_3{x}_3+{\beta }_4{x}_4+ϵ,$$

(1)

where $x_1,x_2,x_3,x_4$ are dummy variables obtained from five levels, y is the response variable and $ϵ$ is a random error generated from a normal distribution with mean 0 and a variance of 1.5. We set ${\beta }_0=2,{\beta }_1=1.5,{\beta }_2=-1.2,{\beta }_3=0.7,\mathrm{and}{\beta }_4=2.0$, respectively. Then the response y can be simulated based on the model (1). We fit a Bayesian linear regression model to the simulated data with y and the four dummy indicators ([5]). For each estimated regression coefficient, we collect the posterior samples from MCMC to draw the corresponding posterior distributions, as shown in Figure 1. In addition to the plots of posterior distribution, we also denote the posterior median and associated 95% credible intervals. These intervals represent the range within which the true value of each regression coefficient lies with 95% posterior probability, given the model and observed data.

The posterior inference is given based on the full posterior distribution of each coefficient. Rather than relying on point estimates and p-values as in frequentist inference, the Bayesian framework allows us to summarize the uncertainty in each parameter using credible intervals and posterior densities. Specifically, we computed the posterior median and 95% credible interval for each coefficient, which are visualized using density plots that highlight both the posterior distribution and the location of central estimates and intervals. The posterior median for ${\beta }_1$, associated with predictor $x_1$, corresponding to the difference between group 2 and the reference group (group 1), is 1.38, with a 95% credible interval of (0.80, 1.98), indicating strong evidence for a positive effect on the response variable. This implies that, the average effect of group 2 is 1.38 times higher than that of group 1. The 95% credible interval supports a substantial positive association and the significance of $x_1$ in the final model. The coefficient ${\beta }_2$ represents the contrast between group 3 and group 1. The posterior median is –0.67, with a 95% credible interval of (–1.27, –0.05). Since the entire interval falls below zero, this provides strong evidence of a negative association, indicating that the average effect of the individuals in group 3 tends to have significantly lower outcomes than those in group 1. $x_2$ is significant as the credible interval doesn’t cover zero. The coefficient ${\beta }_3$ corresponds to the difference between group 4 and group 1. The posterior median was 0.92, with a 95% credible interval of (0.34, 1.55), indicating that the average effect in group 4 is higher than in group 1. The coefficient ${\beta }_4$ captures the contrast between group 5 and group 1. The 95% credible interval is (1.75, 2.92), indicating a strong and statistically significant positive association. The entire credible interval lies well above zero, which suggests that individuals in group 5 have, on average, substantially higher outcomes than those in group 1. The width of the interval also suggests that the magnitude of this effect is relatively stable across the posterior distribution.

To examine the impact of different coding strategies, we fit Bayesian linear regression models using each of the four coding schemes, enabling a direct comparison of their effects on parameter estimates and interpretations. The model coefficients are given in Table 5.

Using the values of $x_1$-$x_4$ from model 1 and the estimated coefficients from Table 5, we get the predicted score for Group 2 for each coding strategy.

$$\begin{array}{cc}& 2.0982+1.3411\times 1-1.0812\times 0+0.2286\times 0+1.7198\times 0=3.4393\left(\mathit{Dummy}\right)\hfill \\ & 2.5409-0.4512\times 0+0.9021\times 1-1.5095\times 0-0.2159\times 0=3.4430\left(\mathit{Deviation}\right)\hfill \\ & 2.1038+1.3418\times 1-2.4029\times 0+1.2956\times 0+1.4903\times 0=3.4456\left(\mathit{Sequential}\right)\hfill \\ & 2.5411-0.5625\times \left(-{\displaystyle \frac{1}{5}}\right)+1.0504\times \left({\displaystyle \frac{3}{4}}\right)-2.0459\times 0-1.4811\times 0=3.4414\left(\mathit{Helmert}\right)\hfill \end{array}$$

We can find that across the four coding schemes, the predicted scores for group 2 from the Bayesian linear regression model are similar but not identical. This is because in Bayesian linear regression, the posterior estimates combine the likelihood with prior distributions on the regression coefficients. When independent priors are assigned to the coefficients, changing the coding scheme alters the scaling and correlation of the predictors in the design matrix, which in turn influences the estimates. As a result, posterior medians and corresponding predicted scores vary slightly across coding strategies. In contrast, in linear regression, the fitted values are invariant to full-rank linear transformations of the predictors. Different coding schemes for the same categorical variable produce design matrices that are linear transformations of one another, so the estimation yields identical predicted scores even though the individual coefficients and their interpretations differ, which has also been confirmed by [1].

3. Bayesian LASSO Regression

Published literature suggests that LASSO regression has been increasingly adopted in social and behavioral science due to its advantage in handling a large number of predictors compared to standard linear regression ([1,17]). In our study, our focus is on investigating regression analysis with categorical predictors using the Bayesian alternative of LASSO, or Bayesian LASSO ([6]). While Bayesian and frequentist methodologies represent distinct statistical paradigms, they are deeply interconnected. To motivate the use of Bayesian LASSO in our analysis, we begin with a brief overview of frequentist linear regression and LASSO, as well as the corresponding uncertainty quantification procedures.

3.1. Linear Regression, LASSO and Statistical Inference

Consider the following linear regression model with n subjects and p predictors:

$$\mathbf{Y}=\mathbf{X}\beta +ϵ,$$

(2)

where $n\times 1$ vector of $\mathbf{Y}$ denotes the response variable, and $\mathbf{X}$ is an $n\times p$ design matrix with the corresponding regression coefficient vector $\beta$. The random error $ϵ$ follows the multivariate normal distribution $N_n(\mathbf{0},{\sigma }^2{\mathbf{I}}_n)$ with unknown variance parameter ${\sigma }^2$.

In classical linear regression, the ordinary least squares (OLS) method is the most widely used approach for estimating regression coefficients $\beta$. It seeks to find an estimated $\widehat{\beta }$ that minimizes the following prediction error in terms of least square loss:

$${\widehat{\beta }}_{LS}=\underset{\beta }{min}\left|\right|\mathbf{Y}-{\mathbf{X}\beta \left|\right|}_2^2,$$

where $\left|\right|·{\left|\right|}_2$ denotes the L2 norm. Such a least square regression only works when the sample size n is larger than the number of predictors p. In other words, we can only obtain a valid least square estimator $\widehat{\beta }$ in the low-dimensional (i.e., large sample size, low dimensionality) settings. Notably, this does not require any assumptions on the distribution of the error terms. However, to develop statistical inference procedures in order to quantify the uncertainty of $\widehat{\beta }$, we need to assume that model errors follow independent and identically normal distributions with mean zero and constant variance ${\sigma }^2$, which eventually leads to statistical inference procedures in terms of p-values and confidence intervals.

Statistical inference plays an important role in linear regression analysis, especially when the predictors are categorical. For example, as discussed in the previous section, in linear regression with a group of binary indicators obtained through the dummy coding, the regression coefficients represent the mean difference between the category of interest and the baseline category. Even if the estimated regression coefficient is non-zero, we cannot conclude that the two groups differ at the population level, as the dataset represents only a small sample drawn from the entire population. We can only claim that the difference is statistically significant if the associated p-value is less than 0.05 or the corresponding 95% confidence interval does not include zero. Quantifying uncertainty is crucial for providing scientifically sound conclusions.

In the presence of a large number of features, linear regression is no longer suitable, especially when the goal is to identify a subset of predictors associated with the response variable y. Therefore, LASSO ([3]) has been developed as a penalized least square regression with $L_1$ penalty of the following form:

$${\widehat{\beta }}_{LASSO}=\underset{\beta }{min}\left|\right|\mathbf{Y}-{\mathbf{X}\beta \left|\right|}_2^2+{\lambda \left|\right|\beta \left|\right|}_1,$$

(3)

where $\left|\right|·{\left|\right|}_1$ denotes the L1 norm, and $\lambda$ is a tuning parameter indicating the amount of shrinkage imposed on regression coefficients $\beta$. LASSO can be viewed as a regularized least squares regression with a constraint on the magnitude of $\beta$ in terms of its L1 norm ${\left|\right|\beta \left|\right|}_1$. When $\lambda$=0, there is no shrinkage on estimating $\beta$, and LASSO reduces to least squares regression. As $\lambda \to \infty$, the constraint on ${\parallel \beta \parallel }_1$ becomes increasingly stringent, leading to more zero-valued components in $\beta$. In other words, the model becomes increasingly sparse. When the regression coefficient, say ${\widehat{\beta }}_j$, is zero, the corresponding jth predictor is no longer associated with y. Therefore, as $\lambda$ increases, model complexity decreases, and fewer predictors are selected in the final model. Choosing an appropriate tuning parameter is thus crucial in LASSO to retain a meaningful subset of predictors ([18]).

Reference [1] have conducted a detailed analysis on the impact of coding strategy choice in LASSO regression with categorical predictors; however, statistical inference procedures of LASSO have not been investigated. In fact, although methods for quantifying uncertainty in LASSO regression have been extensively developed in the literature, such as in [19,20,21,22], among many others ([7]), they are primarily grounded in theoretical studies, making them difficult for practicing statisticians and researchers in the behavioral and social sciences to apply in substantive research. This discrepancy between the availability of frequentist LASSO inference procedures and their limited application in the behavioral and social sciences has motivated us to explore the Bayesian alternative to LASSO in this study.

3.2. From LASSO to Bayesian LASSO

We first illustrate the major difference between frequentist and Bayesian realms using the linear regression model outlined in (2). In the frequentist framework, the model parameters of model (2), consisting of regression coefficient $\beta$ and variance parameter ${\sigma }^2$, are treated as fixed but unknown constants. For $\beta$, we fit a least square regression model to obtain its estimate ${\widehat{\beta }}_{LS}$ and then derive the inference measures in terms of p-value and confidence intervals for uncertainty quantification. On the other hand, within the Bayesian framework, all model parameters, including $\beta$, are treated as random variables, with a prior distribution placed on them. We then follow the Bayes’ theorem to derive its entire posterior distributions ([5]):

$$\mathrm{Posterior}\mathrm{Distribution}\left(\beta \right)\propto \mathrm{Likelihood}(\beta ;\mathbf{Y},\mathbf{X})\times \mathrm{Prior}\left(\beta \right),$$

(4)

where the full posterior distribution of $\beta$ can be obtained via posterior sampling using Markov Chain Monte Carlo (MCMC). This fully Bayesian analysis enables us not only to derive point estimates, such as the posterior mean, median, or any percentile of interest, but also to conduct inference through Bayesian credible intervals or false discovery rates (FDR).

When applying Bayesian alternatives of LASSO to deal with a large number of predictors is of interest, it is crucial to develop an appropriate shrinkage prior on $\beta$ to induce Bayesian LASSO. The following independent and identical Laplace prior has been proposed on ${\beta }_j$, ($j=1,\dots ,p$) ([6]),

$$\pi \left({\beta }_j\right|{\sigma }^2)={\displaystyle \frac{\lambda }{2\sqrt{{\sigma }^2}}}{e}^{-\lambda |{\beta }_j|/\sqrt{{\sigma }^2}},$$

(5)

which is conditional on ${\sigma }^2$ and thus leads to the posterior distribution with unimodality. With the normality assumption on model errors that $ϵ\sim N_n(\mathbf{0},{\sigma }^2{\mathbf{I}}_n)$, the likelihood function of model (2) can be concisely expressed as :

$$f\left(\mathbf{Y}\right|\beta ,{\sigma }^2)\propto \exp \{-{\displaystyle \frac{1}{2{\sigma }^2}}\left|\right|\mathbf{Y}-{\mathbf{X}\beta \left|\right|}^2\},$$

(6)

where ∝ denotes that the likelihood function $f\left(\mathbf{Y}\right|\beta ,{\sigma }^2)$ is proportional to the exponential kernel up to certain normalizing constant not involving $\beta$. We can then derive the posterior distribution of $\beta$ following the Bayes rule (4) by multiplying the conditional Laplace prior (5) from 1 to p across j to the likelihood function (6),

$$\pi (\beta ,{\sigma }^2|y)\propto \exp \{-{\displaystyle \frac{1}{2{\sigma }^2}}\left|\right|\mathbf{Y}-{\mathbf{X}\beta \left|\right|}_2^2-\lambda \sum _{j=1}^p|{\beta }_j\left|\right\}.$$

(7)

The connection between LASSO and its Bayesian counterpart can be clearly revealed by comparing the penalized least square formulation in (3) and its posterior distribution in (7). Minimizing the penalized least square loss is equivalent to maximizing (7). Therefore, the frequentist LASSO estimate is equivalent to the corresponding maximum a posteriori (MAP) estimate (i.e., Bayesian posterior mode estimate) under the Bayesian framework. [6] have shown that the conditional Laplace prior on ${\beta }_j$ defined in (5) guarantees that the resulting posterior mode is unique.

As we discussed in Section 3.1, the theoretical nature of inference with LASSO has made it less accessible to practitioners in the social and behavioral sciences. In contrast, fully Bayesian analysis leverages posterior samples drawn from MCMC to conduct posterior inference on model parameters. As long as practitioners are familiar with standard Bayesian analysis, they can run Bayesian LASSO and readily use the posterior samples to conduct inference on $\beta$. LASSO regression with categorical predictors has been carefully examined in [1]. However, the lack of conducting statistical inference with LASSO, due to the aforementioned theoretical challenges, has motivated us to explore bridging this gap using the Bayesian LASSO.

4. Real Data Analysis

To explore the potential impacts of coding strategy on important predictors, we analyze the 5-year real-world MEPS data using Bayesian LASSO, LASSO and linear regression. The MEPS data set has a sample of 98,163 participants. In the analysis, the response of interest is the PCS score. We consider six categorical predictors: MS (2 categories), Sex (2 categories), Race (4 categories), Region (4 categories), Education (5 categories), Insurance (3 categories), and two continuous variables: Age and ECI (Number of Elixhauser comorbidity conditions). The performance of three models under comparison has been examined under different coding strategies in terms of variable selection, prediction, and inference procedures.

4.1. Variable Selection

We investigate whether the choice of coding strategy influences the variable selection with Bayesian LASSO. The results are provided in Table 6, which shows that selected variables vary under different coding strategies. For example, variable Race corresponds to three predictors after conversion. With dummy coding, all three race binary indicators are included in the model. With deviation coding, the predictor measuring the difference between non-Hispanic black and the average is excluded from the final model. However, the sequential coding strategy excludes the predictor indicating the difference between Hispanic and other, while the Helmert coding strategy excludes the predictor representing the difference between non-Hispanic black and other. From this example, we can find that a researcher may conclude that the PCS scores differ across the race of Hispanic and other with the dummy coding strategy, while using the sequential coding strategy, the opposite conclusions will be made. Variable selection with different coding strategies has also been assessed using LASSO and linear regression, as shown in Table A3 and Table A4 in the Appendix, respectively. Notice that linear regression alone cannot perform variable selection. For a direct comparison to the other two methods. we compute the 95% confidence interval for each regression coefficient, and exclude predictors if the corresponding confidence intervals cover zero. We can find that the important predictors identified through linear regression and LASSO vary depending on the coding strategy, as different coding schemes alter the model’s parameterization and thus the estimated regression coefficients.

A cross-comparison of results from Table 6, Table A3 and Table A4 indicates that on MEPS data, Bayesian LASSO leads to exclusion of variables under all four types of coding. The number of excluded features under Dummy, Deviation, Sequential, and Helmert coding schemes are 2, 1, 1, and 2, respectively. In comparison, linear regression excludes 0, 1, 1, and 2 features under the same codings. The selection results between the two methods are more similar, as LASSO eliminates 0, 1, 0, and 1 variables, respectively. Although it may initially seem surprising that the variable selection results from Bayesian LASSO and linear regression are more similar to each other, a closer examination provides a clear explanation. Both Bayesian LASSO and linear regression rely on inference-based measures—such as confidence or credible intervals—for feature selection. In contrast, LASSO selects variables by applying regularization and eliminating features whose coefficient estimates are exactly zero. Overall, it can be concluded that the variable selection in both Bayesian LASSO, linear regression and LASSO is influenced by the choice of coding strategies for categorical variables.

4.2. Prediction Accuracy

Under different coding strategies, we have assessed the prediction performance for all methods under comparison in terms of (1) predicted category scores and (2) the least absolute deviation (LAD) error, which is defined as $LAD={\displaystyle \frac{1}{n}}{\sum }_{j=1}^p|Y_i-\widehat{Y_i}|$. Choosing Education as an example, we examine whether the predicted score for each education level is the same under different coding strategies using Bayesian LASSO. Table 7 shows predicted PCS scores using Bayesian LASSO corresponding to the five levels of Education, with the last column representing the actual category means. It can be observed that under different coding strategies, the predicted scores are of the same magnitude for the same category, with slight differences. Although the predicted scores are close to the actual category means, there are no models where the predicted scores are equal to the true category means. Similar patterns can be found with LASSO in Table A2 in the Appendix. However, the predicted category scores are exactly the same when using linear regression, which are also equal to the actual category means, as shown in Table A1 in the Appendix. Mathematically, any coding transformation corresponds to a change in the coefficient estimates but maintains the same fitted values in linear regression, as the full model space is explored without restriction.

Next, we assess predictive performance by computing the previously defined LAD error for all the three models by including all categorical and continuous variables under different coding strategies. The results are listed in Table 8. In the case of linear regression, the LAD error remains unchanged across different coding strategies because the predicted values are invariant to how categorical variables are encoded. However, for Bayesian LASSO, shrinkage effects are not the same under different coding schemes, leading to different penalized estimates and, consequently, varying prediction errors in terms of LAD. Under the 4 parameterizations, the model performance is impacted under Bayesian LASSO and LASSO, while remaining the same with linear regression.

4.3. Statistical Inference

Inference procedures play a crucial role in statistical analysis. However, frequentist high-dimensional methods, including LASSO, typically rely on complex asymptotic theory to develop inference procedures ([23]). This reliance poses significant challenges for practitioners seeking to understand and apply these methods to real-world data analysis. A more detailed discussion of the obstacles associated with implementing frequentist inference procedures in practical settings can be found in [7]. On the other hand, Bayesian approaches overcome this difficulty by providing a principled and coherent framework of conducting statistical inference through posterior sampling regarding all model parameters ([5]). By building up Bayesian hierarchical models that leverage strength from prior information and observed data, fully Bayesian analysis can characterize the entire posterior distribution of model parameters via sampling based on Markov Chain Monte Carlo (MCMC) and techniques alike. Therefore, uncertainty quantification measures including standard summary statistics (such as median, mean and variance), credible intervals and posterior estimates with false discovery rates (FDR) can be readily obtained.

We perform statistical inference in terms of marginal Bayesian credible intervals with Bayesian LASSO on the MEPS dataset under the four coding systems. Figure 2 shows posterior distributions of three regression coefficients resulting from converting Race following dummy coding strategy. Using Hispanic individuals as the reference group, for Non-Hispanic white individuals, the posterior median is -0.45 with a 95% credible interval ranging from -0.58 to -0.31. Since this interval does not include zero, it suggests a statistically significant negative association compared to Hispanic individuals, indicating that Non-Hispanic whites have a significantly lower level of PCS score than the reference group. For Non-Hispanic black individuals, the posterior median is -0.73 with a 95% credible interval of -0.90 to -0.53. This interval also excludes zero, providing strong evidence of a significant negative difference in the outcome relative to the Hispanic group. Similarly, for individuals classified as other, the posterior median is -0.68, and the 95% credible interval extends from -0.90 to -0.46. This interval likewise does not include zero, indicating a statistically significant negative association with PCS score compared to Hispanic individuals. Taken together, these results demonstrate that all three racial groups exhibit significantly lower PCS scores compared to the Hispanic reference group. As the posterior credible intervals exclude zero, it suggests these predictors should be kept in the final model, which corresponds to the variable selection results in Table 6.

The posterior median, 95% credible interval and the corresponding interval length of regression coefficients associated with all predictors under the four coding strategies using Bayesian LASSO have been shown in Table A5 in the Appendix. The covariates with the credible interval excluding zero are considered to be significantly associated with the outcome and are included in the final model, which is listed in Table 6. We can find that while the posterior medians and credible intervals of regression coefficients vary across different coding strategies, their overall magnitudes are close. Although the underlying relationships between the predictors and outcome remain the same across coding strategies, the parameterization changes. Thus, while the parameter estimates and their associated uncertainty intervals differ, the substantive conclusions about the strength and direction of associations remain relatively stable. This explains why the posterior medians and credible intervals vary numerically but reflect comparable effect sizes across coding strategies. This suggests that although the numerical representations of the categorical variables influence the specific parameter estimates and their associated uncertainty ranges, the underlying effect sizes are consistently captured across coding strategies. These differences in posterior summaries highlight the sensitivity of Bayesian LASSO inference to the choice of coding strategies, yet the comparable magnitude of estimates indicates that the substantive conclusions about variable importance remain relatively stable. As a comparison, we present estimation results along with information on confidence intervals by fitting linear regression in Table A6 in the Appendix. Table A7 in the Appendix shows the estimation results using LASSO, which is similar to the analysis conducted in [1] where only point estimates of LASSO have been considered. In the linear regression analysis, the estimated coefficients remain consistent across different coding strategies. In contrast, the LASSO approach exhibits greater sensitivity to coding choices and certain coefficients are shrunk exactly to zero under some strategies, while the same coefficients are retained as nonzero under others. Unlike Bayesian methods, both linear regression and LASSO provide only point estimates of the coefficients, without a direct framework for quantifying uncertainty.

4.4. Convergence

Assessing convergence of Markov Chain Monte Carlo (MCMC) is critical in Bayesian analysis, as the validity of all Bayesian posterior estimates and inferences depends on that MCMC has converged and reached the corresponding stationary distribution ([5]). If the chains do not converge properly, the posterior samples drawn from may not accurately characterize correct underlying posterior distributions, leading to biased estimates and inferences, as well as misleading conclusions ([24]). To assess the reliability of posterior estimates, we evaluated standard MCMC convergence diagnostics for all model parameters with the potential scale reduction factor (PSRF) ([25]). By running multiple MCMC chains on the same dataset, PSRF compares the variance within chains to the variance between chains. If all chains have converged to the target posterior distribution, these variances should be similar, leading to a value close to 1. Otherwise, the chains do not mix well if PSRF value is much larger than 1.

In this study, we use PSRF ≤ 1.1 ([26]) as the cut-off point which indicates that chains converge to a stationary distribution. In MCMC, the initial iterations are often affected by the choice of starting values and may not adequately represent draws from the stationary posterior distribution. To mitigate this influence and reduce bias in posterior estimation, these early samples are commonly discarded as burn-in iterations ([5]). In this study, the Gibbs sampling is implemented with 10,000 iterations with the first 5000 as burn-ins. The convergence of the MCMC chains after burn-ins has been checked for all predictors under the four coding strategies, which can be found in Figure A1, Figure A2, Figure A3 and Figure A4 in the Appendix. Take Race from dummy coding as an example in Figure 3. The PSRF trajectories for the three Race dummy variables corresponding to Non-Hispanic white, Non-Hispanic black, and other race categories, respectively. For each parameter, the median PSRF (solid black line) rapidly approaches and stabilizes near 1.00, and the upper 97.5% quantile (dashed red line) remains well below the commonly used threshold of 1.1 after early iterations. This pattern indicates good mixing and convergence of the Markov chains for all three race coefficients. The initial variability observed in the early iterations diminishes quickly, further supporting the conclusion that the chains have converged. These results suggest that posterior inference for the race-related coefficients is reliable and based on well-converged MCMC samples.

5. Discussion

In this study, we investigate the impact of coding strategies on categorical predictors in variable selection, prediction, and statistical inference with Bayesian LASSO in MEPS data. Comparisons against frequentist approaches such as linear regression and LASSO have also been performed. Our study complements existing work, such as [1], by adopting a Bayesian framework and incorporating formal statistical inference procedures into the analysis. Huang et al. ([1]) have also evaluated group LASSO as an alternative to LASSO on regression with categorical predictors, and concluded that group LASSO leads to overfitting few instead of all dummy predictors within the same group are needed. While we agree with their conclusion regarding the application of group LASSO, we want to point out that sparse group LASSO ([27,28]) is a promising regularization method that seeks to achieve sparsity on both group level and within group level. Therefore, sparse group LASSO type of methods are potentially promising in the scenario when selection of important categorical predictors within groups are of interest.

Due to heterogeneity of complex diseases such as multiple sclerosis and cancers, disease phenotypes of interest usually follow heavy-tailed distributions and have outlying observation. Therefore, robust statistical methods, especially robust regularization and variable selection methods that can safeguard against outliers and skewed distributions are demanded ([18]). Recently, the advantages of robust Bayesian variable selection methods over their frequentist counterparts, particularly in the context of statistical inference, have been investigated in [7]. It will be interesting to explore how robust Bayesian analysis can facilitate modeling with categorical predictors. For example, the robust Bayesian sparse group LASSO model proposed in [29] offers uncertainty quantification, which is typically unavailable in corresponding frequentist approaches ([27,28]). In future work, we also plan to extend the current methodology to other types of phenotypic traits, such as survival and longitudinal outcomes.