Discussion: Bayesian regularization and effect smoothing for categorical predictors

Abstract

Keywords

bayesian regularization spike and slab Prior bayesian effect fusion MCMC

1 Introduction

We would first like to thank Gerhard Tutz and Jan Gertheiss for their profound review of regularization methods for categorical variables, either as covariates or as response variables in regression models. Categorical variables are rather the rule than the exception in regression analyses, particularly in the medical, social and economic sciences. It is amazing that it took quite a long time from ridge regression (Hoerl and Kennard, 1970) to versions of regularization methods, which take into account the specific structure of categorical covariates. As Tutz and Gertheiss predominantly discuss penalized estimation, we will focus on the Bayesian perspective to regularization and effect fusion for categorical covariates.

Regularization and sparsity in regression type models have been addressed from a Bayesian point of view in the literature, see for example, Fahrmeir et al. (2010) for an overview on Bayesian regularization and (George and McCulloch 1995), (Ishwaran and Rao 2005) and (O'Hara and Sillanpäa 2009) for approaches to Bayesian variable selection. However, the specific issues arising for categorical covariates have not yet received much attention, with the exception of Chipman (1996), who considers a Bayesian approach for groupwise selection of level effects of a categorical covariate (the Bayesian pendant to Section 3.2.1 of Tutz and Gertheiss 2016) and the recently proposed Bayesian pendant to smoothing predictors and effect fusion based on spike and slab priors in Pauger and Wagner (2016).

Here, we will discuss the relation between regularization penalties and prior distributions and then describe a Bayesian approach for smoothing of level effects of categorical covariates. Finally, we will briefly describe spike and slab prior distributions which encourage a sparse representation of their effects.

2 Bayesian regularization of effects of categorical covariates

Specification of a Bayesian linear regression model requires not only a model for the data, for example, the linear regression model (3) of Tutz and Gertheiss (2016),

y = α + \sum_{j = 1}^{p} X_{j} β_{j} + ε, ε \sim N (0, σ^{2} I),

(2.1)

but also a prior distribution $p (Θ)$ on the model parameters $Θ$ , that is, $Θ = (α, β, σ^{2})$ . The joint posterior distribution is then given as

p (Θ | y) \propto p (y | Θ) p (Θ),

respectively taking logarithms as

l (Θ | y) = l (y | Θ) + l (Θ) + c .

Obviously, the log-posterior corresponds to the penalized log-likelihood with

J_{λ} (Θ) = - l (Θ) - c .

Consider, for example, a prior distribution of the structure $p (α, β, σ^{2}) \propto p (β)$ with a Normal distribution for the regression coefficients, $β \sim N (0, I τ^{2})$ . Then the log-posterior is given as

l (β | y) = l (y | β) - \frac{1}{2 τ^{2}} | | β | |_{2}^{2}

and the Bayesian posterior mode estimate is equivalent to the regularized ML-estimate with

λ = \frac{1}{2 τ^{2}}

Bayesian regularization is based on specifying a prior distribution $p (β | ξ)$ where the hyperparameters $ξ$ include parameters controlling complexity, shrinkage or smoothing, see Fahrmeir et al. (2010). An advantage of Bayesian regularization is that no cross-validation is needed, as these hyperparameters are estimated jointly with the parameters of interest and, hence, uncertainty on the regularization parameters can also be assessed.

Therefore, to achieve Bayesian regularization for the effects of categorical covariates, appropriate prior distributions for the coefficients $β_{j}$ are needed. As pointed out in Tutz and Gertheiss (2016), to smooth effects of categorical covariates, it is essential to take into account their scale: For ordinal covariates, coefficient smoothing is equivalent to shrinking effect differences of adjacent categories to zero. In contrast, there is no natural order of the categories for a nominal covariate and, therefore, smoothing of coefficients implies to shrink all pairwise effect differences to zero. We discuss prior distributions for Bayesian regularization based on the appropriate effect differences for ordinal covariates in Section 2.1 and for nominal covariates in Section 2.2. To specify these prior distributions, following Tutz and Gertheiss (2016), we assume that the effect of the reference category is defined as $β_{j 0} = 0$ .

2.1 Ordinal covariates

According to Tutz and Gertheiss (2016), Section 3.1.1., if all $p$ covariates are ordinal, smoothing of the effects $β_{jr}$ of adjacent categories $r = 1, \dots, k_{j}$ can be enforced by a quadratic penalty on the effect differences $θ_{j, r} = β_{jr} - β_{j, r - 1}$ ,

J (θ) = \sum_{j = 1}^{p} \sum_{r = 1}^{k_{j}} θ_{j, r}^{2} .

(2.2)

Here, $θ = (θ_{1}, \dots, θ_{p})$ is the vector of all effect differences and $θ_{j} = (θ_{j, 1}, \dots, θ_{j, k_{j}})$ denotes the effect differences of the $j$ -th covariate. The quadratic penalty (2.2) corresponds to the specification of independent Normal priors on all $θ_{j, r}$ which share the scale parameter $τ^{2}$ , i.e. $θ_{j, r} \sim N (0, τ^{2}) .$

To perform Bayesian effect smoothing, it is natural to allow for covariate specific scale parameters. We assume prior independence across covariates $p (θ) = \prod_{j = 1}^{p} p (θ_{j})$ and specify the prior on $θ_{j}$ hierarchically as

\begin{matrix} θ_{j} & \sim N (0, τ_{j}^{2} I_{j}) \end{matrix}

(2.3)

τ_{j}^{2} ~ G^{- 1} (v_{j}, Q_{j}),

(2.4)

where

I_{j}

is a unit matrix of dimension

k_{j}

and

G^{- 1}

denotes the inverse Gamma distribution.

The reparameterized effects $θ_{j}$ are the regression effects for split-coded levels of the $j$ -th covariate. Both $θ_{j}$ and the regression effects $β_{j}$ for dummy-coded levels are of the same dimension and uniquely related by

β_{j} = E_{j} θ_{j},

where

E_{j}

is a lower triangular matrix of dimension

k_{j} \times k_{j}

with elements

1

. Hence, the resulting prior distribution on

β_{j}

N (0, B_{0 j})

, where

B_{0 j} = τ_{j}^{2} E_{j} E_{j}^{'} .

2.2 Nominal covariates

As noted above, for nominal covariates, no natural order of the categories exists and, hence, smoothing of category effects implies shrinking of all effect differences

θ_{j, rs} = β_{jr} - β_{js}

to zero. To construct such a prior, we denote by

θ_{j} = (θ_{j, 10}, \dots, θ_{j, k_{j}, k_{j - 1}})

the vector of all effect contrasts. The first

k_{j}

elements of

θ_{j}

are the effect contrasts with respect to the reference category and correspond to

β_{j}

and hence

θ_{j}^{'} = (β_{j}^{'}, ϑ_{j}^{'})^{'}

, where

ϑ_{j}

is the vector of all effect contrasts not involving the reference category. Elements of

θ_{j}

are linearly dependent, as

θ_{j, rs} = β_{jr} - β_{js} = θ_{j, r 0} - θ_{j, s 0} .

Hence, the smoothing prior for $β_{j}$ can be constructed from a prior distribution on all effect differences $θ_{j}$ by taking into account the linear restrictions

A_{j} β_{j} - ϑ_{j} = 0_{j},

(2.5)

where

A_{j}

is the appropriate

\frac{k_{j} (k_{j} - 1)}{2} \times k_{j}

coefficient matrix and

0_{j}

is the corresponding vector of zeros.

Slightly deviating from the corresponding prior for an ordinal covariate, the prior on the effect differences is specified as

\begin{matrix} θ_{j} & \sim N (0, γ_{j} τ_{j}^{2} I), \\ τ_{j}^{2} & \sim G^{- 1} (ν_{j}, Q_{j}), \end{matrix}

where

γ_{j} = \frac{k_{j} + 1}{2}

takes into account the number of levels. As shown in Pauger and Wagner (2016), the resulting prior on

β_{j}

under the linear restrictions (2.4) is

N (0, B_{0 j})

where

B_{0 j} = γ_{j} τ_{j}^{2} (I_{j} - A_{j}^{'} (A_{j} A_{j}^{'} + I)^{- 1} A_{j}) .

2.3 Posterior inference

In a linear regression model with different types of covariates (metric, ordinal, nominal) assuming prior independence across covariates, the resulting prior distribution for the vector of all regression coefficients $β$ is a Normal distribution $N (0, B_{0})$ , where $B_{0}$ is block-diagonal with elements $B_{0 j}$ . With standard priors for the intercept $α \sim N (0, M)$ and the error variance $σ^{2} \sim G^{- 1} (s_{0}, S_{0})$ , posterior inference using Markov chain Monte Carlo (MCMC) methods is straightforward and requires only standard Gibbs sampling steps:

Sample $β$ from the full conditional multivariate Normal distribution.

Sample the error variance $σ^{2}$ from the full conditional inverse Gamma distribution.

For $j = 1, . . ., p$ , sample the scale parameter $τ_{j}^{2}$ from its full conditional inverse Gamma distribution.

2.4 Application

To illustrate Bayesian effect smoothing, we reanalyze the Munich rent standard data used in Gertheiss and Tutz (2010). The response variable monthly rent per square meter is modelled as a function of five binary (V1–V5), four ordinal (V6: ‘year of construction’, V7: ‘number of rooms’, V8: ‘quality of residential area’ and V9: ‘floor space’ in ordered classes) and one nominal covariate (V10: ’urban district’ with 25 categories), see Gertheiss and Tutz (2010) for a more detailed description of the data set. To perform Bayesian analysis, we choose priors for the regression coefficients as described in Sections 2.1 and 2.2 and uninformative priors for the intercept $α \sim N (0, 10 000)$ and the error variance $σ^{2} \sim G^{- 1} (0, 0)$ . We set $ν_{j} = 5$ and vary the scale parameters $Q_{j}$ in a wide range. Figure 1 shows paths of the resulting coefficients $β_{j}$ as a function of the corresponding scale parameter $Q_{j}$ .

Figure 1:

Coefficient paths for the ordinal predictors ’year of construction’, ’number of rooms’, ’quality of residential area’ and ’floor space’ and the nominal predictor ’urban district’.

As the expected prior variance of the effect differences decreases with $Q_{j}$ , we see increasing shrinkage of effect differences $θ_{j, rs}$ to zero and smoothing of the level effects. However, different from smoothing and fusion penalties, except for very small values of $Q_{j}$ , effect differences are not exactly zero.

3 Spike and slab priors

Sparsity in the representation of the effects of a categorical covariate can be achieved by replacing the Normal prior on the elements of $θ_{j}$ by spike and slab priors. These priors are a mixture of two components: a spike at zero which shrinks small coefficients aggressively to zero and a slab component with its mass spread over a range of plausible values, which induces only minor shrinkage. Spike and slab priors are specified hierarchically by introducing indicator variables $δ_{j, rs}$ which take the value one, if $θ_{j, rs}$ is assigned to the slab component. The prior on $θ_{j, rs}$ is given as

\begin{matrix} p (θ_{j, rs} | δ_{j, rs}, ω_{j}, τ_{j}^{2}, γ_{j}) & \sim δ_{j, rs} N (0, τ_{j}^{2} γ_{j}) + (1 - δ_{j, rs}) N (0, \frac{τ_{j}^{2}}{r} γ_{j}), \end{matrix}

(3.1)

p (δ_{j, r s} = 1) = ω_{j},

(3.2)

where

r

is a large value. The inclusion probability

ω_{j}

can be fixed or assigned a Beta prior

ω_{j} \sim B (a_{j 0}, b_{j 0})

. We subsume all indicators associated with the

j

-th covariate in the vector

δ_{j}

For a nominal covariate, we derive in Pauger and Wagner (2016) the structure of the prior variance matrix $B_{0 j} (δ_{j})$ of $β_{j}$ and show that $δ_{j, rs} = 0$ implies that the prior partial correlation of $β_{jr}$ and $β_{js}$ is close to $1$ .

Posterior sampling requires to add only the following two simple sampling steps to the scheme described in Section 2.3, if a hyperprior is specified for $ω_{j}$ :

For $j = 1, \dots, p,$ sample the binary indicators $δ_{j, rs}$ from $P (δ_{j, rs} = 1 | y, β_{j}, τ_{j}^{2}, ω_{j})$ and update the prior variance matrix $B_{0 j} (δ_{j}) .$

For $j = 1, \dots, p,$ sample $ω_{j}$ from its full conditional Beta distribution;

otherwise sampling

ω_{j}

is omitted.

To perform Bayesian effect fusion for the Munich rent data, we set $r = 30 000$ , $ω = 0.5$ , $ν = 5$ and $Q = 2$ for all covariates. We run MCMC for 20 000 iterations after a burnin of $10 000$ . Level effects are fused based on an estimated posterior inclusion probability ${\hat{p}}_{j, rs} > 0.5$ and the resulting model (Model 1) is reported in Table 1.

Table 1:

Munich rent data, Model 1: Model selected by Bayesian effect fusion (grouping of level effects is indicated by spaces)

V6	0 1 2 3 4 5 6 7 8 9
V7	0 1 2 3 4 5
V8	0 1 2
V9	0 1 2 3 4 5 6 7 8 9 10 11 12
V10	0 1 2 3 4 8 11 12 17 5 6 7 9 10 13 14 15 16 18 19 20 21 22 23 24

For covariates $V 7$ and $V 8$ , the same levels as in the model selected in Gertheiss and Tutz (2010), which we will denote by Model (2) are fused. However, fusion is more pronounced in the Bayesian version for the ordinal covariates $V 6$ (7 vs. 8 groups) and $V 9$ (6 vs. 7 groups), and particularly for the nominal covariate $V 10$ , with only two groups under Bayesian fusion and 10 groups using regularization, see also Table 2.

Table 2:

Munich rent data, Model 2: Model selected by regularized regression (grouping of level effects is indicated by spaces)

V6	0 1 2 3 4 5 6 7 8 9
V7	0 1 2 3 4 5
V8	0 1 2
V9	0 1 2 3 4 5 6 7 8 9 10 11 12
V10	0 2 1 3 4 11 17 12 8 5 7 9 14 16 18 19 20 24 6 10 22 13 15 21 23

Finally, in Table 3, we compare the full model, that is, the model with separate regression effects for each covariate level to refits of the two models selected under the different approaches. Both Models 1 and 2 have considerably fewer parameters as well as lower values of the information criteria Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC) than the full model. AIC is smaller for the model obtained by regularization (Model 2), whereas BIC is smaller for the model selected by the Bayesian approach (Model 1), which illustrates that the Bayesian version of effect fusion implies a higher penalty.

Table 3:

Munich rent data: Model selection criteria for different models

	Number of coefficients	AIC	BIC
Full Model	58	8 651.4	8 983.4
Model 1 (Bayesian Fusion)	22	8 616.5	8 745.9
Model 2 (Fusion Penalty)	32	8 604.3	8 789.9

Footnotes

Acknowledgments

This research was supported by the Austrian Science Fund (FWF), P25850.

References

Chipman

(1996) Bayesian variable selection with related predictors. The Canadian Journal of Statistics , 24, 17–36.

Fahrmeir

Kneib

Konrath

(2010) Bayesian regularisation in structured additive regression: A unifying perspective on shrinkage, smoothing and predictor selection. Statistics and Computing , 20, 203–19.

George

McCulloch

(1995) Approaches for Bayesian variable selection. Statistica Sinica , 7, 339–73.

Gertheiss

Tutz

(2010) Sparse modelling of categorical explanatory variables. The Annals of Applied Statistics , 4 (4), 2150–80.

Hoerl

Arthur E

Kennard

Robert W

(1970) Ridge regression: Biased estimation for nonorthogonal problems. Technometrics , 12, 55–67.

Ishwaran

Rao

(2009) Spike and slab variable selection: Frequentist and Bayesian strategies. Annals of Statistics , 33, 730–73.

O'Hara

Sillanpäa

(2009) A review of Bayesian variable selection methods: What, how and which. Bayesian Analysis , 4, 85–117.

Pauger

Wagner

(2016) Bayesian Effect Fusion for Categorical Predictors. Technical Report, Linz: Johannes Kepler University.

Tutz

Gertheiss

(2016) Regularized regression for categorical predictors. Statistical Modelling , in press.