# Bayesian Modal Regression based on Mixture Distributions

Qingyang Liu
Department of Statistics
University of South Carolina
Columbia, SC 29201
<EMAIL_ADDRESS>
& Xianzheng Huang
Department of Statistics
University of South Carolina
Columbia, SC 29201
<EMAIL_ADDRESS>
& Ray Bai
Department of Statistics
University of South Carolina
Columbia, SC 29201
<EMAIL_ADDRESS>

###### Abstract

Compared to mean regression and quantile regression, the literature on modal
regression is very sparse. We propose a unified framework for Bayesian modal
regression based on a family of unimodal distributions indexed by the mode
along with other parameters that allow for flexible shapes and tail behaviors.
Following prior elicitation, we carry out regression analysis of simulated
data and datasets from several real-life applications. Besides drawing
inference for covariate effects that are easy to interpret, we consider
prediction and model selection under the proposed Bayesian modal regression
framework. Evidence from these analyses suggest that the proposed inference
procedures are very robust to outliers, enabling one to discover interesting
covariate effects missed by mean or median regression, and to construct much
tighter prediction intervals than those from mean or median regression.
Computer programs for implementing the proposed Bayesian modal regression are
available at https://github.com/rh8liuqy/Bayesian_modal_regression.

_Keywords_ Mode $\cdot$ Heavy-tailed distribution $\cdot$ Outlier $\cdot$
Unimodal distribution

##  1 Introduction

There is an abundance of literature on mean regression models which model the
conditional mean of a response variable $Y$ given a set of covariates
$\boldsymbol{X}$. However, it is no secret that the mean is sensitive to
outliers. Median regression – or more generally, quantile regression – is
robust to outliers and is thus an appealing alternative to mean regression
(Koenker et al.,, 2017). Besides the mean and median, the mode is yet another
commonly used measure of central tendency. Compared with mean or median
regression, modal regression concerns the conditional mode of $Y$ given
$\boldsymbol{X}$ and is much less explored (Sager and Thisted,, 1982; Lee,,
1989, 1993), especially in the parametric framework.

But why are modal regression models useful additions to the well-established
mean and median regression models? For unimodal and asymmetric distributions,
intervals around the conditional mode typically have higher coverage
probability than intervals of the same length around the conditional mean or
median (Yao and Li,, 2014; Xiang and Yao,, 2022). Consequently, prediction
intervals from modal regression tend to be narrower than those for mean or
median regression when data arise from a unimodal and skewed distribution. By
construction, modal regression explores the relationship between the “most
probable” value of $Y$ given $\boldsymbol{X}$, and thus offers a highly
interpretable representative value of the response. Thanks to the nature of
the mode, modal regression is extremely robust to outliers that can obscure
some inherent covariate effect suggested by the majority of observations,
making it a worthy rival of median regression as an alternative to mean
regression in regard to feature discovery.

A major challenge in building parametric modal regression models is
constructing an appropriate distribution family that subsumes asymmetric,
symmetric, light-tailed, and fat-tailed distributions. In this paper, we
propose the general unimodal distribution (GUD) family, which is a subfamily
of the general two-component mixture distribution family described in Section
3. Members of the GUD family have a location parameter as the mode, in
addition to shape and scale parameters that control the skewness and tail
behaviors. Thus, our framework is appropriate for both asymmetric and
symmetric conditional distributions, as well as both light-tailed and heavy-
tailed distributions. In the extreme case, our framework can also model data
from distributions without any finite moments, which we introduce in Section
3.3. We propose to estimate the conditional mode and the shape/scale
parameters using a Bayesian approach. By placing appropriate prior
distributions on model parameters, our modal regression models can be
implemented straightforwardly using Markov chain Monte Carlo (MCMC) and
provide natural uncertainty quantification through the posterior
distributions.

###  1.1 Existing Work on Modal Regression

Frequentist nonparametric modal regression has been the mainstream in the
limited existing literature on modal regression (see Chen, (2018) for a
comprehensive review). The higher statistical efficiency and greater
interpretability of covariate effects under a parametric framework motivate
some recent development in frequentist parametric modal regression. For
example, Aristodemou, (2014) and Bourguignon et al., (2020) proposed a
parametric modal regression model based on a gamma distribution for a positive
response; Zhou and Huang, (2020) proposed two parametric modal regression
models for a bounded response. Menezes et al., (2021) give a nice review on
these and other parametric modal regression models for a bounded response. In
contrast to these existing parametric modal regression models for bounded
data, the modal regression models in the present manuscript are based on a new
GUD family whose support is the entire real line. Furthermore, our work deals
with Bayesian inference for modal regression.

The literature on Bayesian modal regression is even more sparse. Yu and
Aristodemou, (2012) proposed a nonparametric Bayesian modal regression model
using Dirichlet process mixtures of uniform distributions. Zhou and Huang,
(2022) proposed a parametric Bayesian modal regression model based on a four-
parameter beta distribution whose support is bounded yet unknown. Damien et
al., (2017) introduced a more flexible parametric form of Bayesian modal
regression using mixtures of triangular densities for a response with an
unknown bounded support. Remaining in the parametric framework, a major
strength of our proposed GUD family is that it naturally facilitates data-
driven learning of the skewness and tails of the underlying distribution
supported on the entire real line, while signifying the mode as the central
tendency measure of the response.

###  1.2 Our Contributions

This paper aims to widen the scope of Bayesian modal regression models and
highlight the advantages of these models through analyses of datasets from
real-life applications in several disciplines. We provide a unified framework
of Bayesian modal regression models based on the GUD family that contains a
large variety of unimodal distributions. We adopt the fully Bayesian approach
via MCMC, so that inference of the conditional mode does not rely on
asymptotic approximations. Our method is shown to provide reliable inference
in small sample sizes. The fully Bayesian approach also comes with a
convenient way of constructing prediction intervals from the posterior
predictive distribution that is approximated using a simple random number
generation algorithm for the GUD family. Indeed, the convenience in data
generation via a data augmentation trick (to be discussed on Section 3) is yet
another advantage of the GUD family. Finally, we exploit the model criterion
known as the Bayesian leave-one-out expected log predictive density for model
selection to help practitioners choose an appropriate distribution for their
final analysis.

The main contributions of this paper can be summarized as follows:

  1. 1.

We propose the GUD family that is suitable for Bayesian modal regression. The
GUD family contains distributions that are symmetric or asymmetric,
(non)normal, and/or fat-tailed.

  2. 2.

We formulate rules of prior elicitation for the GUD family. In particular, we
place a flat prior on regression coefficients, weakly informative priors on
all other model parameters, and establish posterior propriety under mild
conditions.

  3. 3.

We provide strategies for constructing prediction intervals and for selecting
an appropriate likelihood for Bayesian modal regression analysis.

  4. 4.

We illustrate the following benefits of our proposed Bayesian modal regression
framework through simulation studies and data applications in economics,
criminology, environmental science, and molecular biology: a) robustness to
outliers, b) more precise prediction, and c) high interpratability of
covariate effects.

The structure of this paper is as follows. In Section 2, we motivate our
proposed Bayesian modal regression framework with two data applications. In
Section 3, we formally define the GUD family and zoom in on several important
members in the family. Section 4 introduces Bayesian inference for these modal
regression models, including prior elicitation, posterior propriety,
uncertainty quantification, and model selection. Section 5 provides simulation
studies that illustrate the strengths of our methodology. In Section 6, we
provide two additional data applications from environmental science and
molecular biology. Section 7 concludes the paper with some remarks about our
Bayesian modal regression framework and some directions for future research.

##  2 Motivating applications

As a prelude to introducing our Bayesian modal regression framework, we first
present results from applying the proposed methodology (to be elaborated in
Sections 3 and 4) to datasets from the economics and the criminology
literature.

###  2.1 Modeling highly right-skewed bank deposits

It is common knowledge to economists that wealth distributions are highly
skewed to the right (Benhabib and Bisin,, 2018). The cumulative nature of
wealth not only has impact on individuals’ net worth, but also has an
influence on assets of large companies, including bank holding companies. In
this example, we analyzed the deposits data of 50 banks and savings
institutions in the United States on July 2, 2010 (Table 3.4.1 in Siegel,
(2016)).

Figure 1 presents the estimated density plot that results from fitting an
intercept-only regression model based on the Double Two-Piece-Student-$t$, or
DTP-Student-$t$, distribution (to be introduced in Section 3) to the dataset,
along with the histogram of the underlying data. From this figure, we can see
that the estimated mode using the DTP-Student-$t$ distribution is close to the
nonparametric mode estimate based on the histogram. This similarity and the
close resemblance of the fitted density to the shape of the histogram indicate
that the DTP-Student-$t$ distribution is an adequate choice for the bank
deposits data.

The other two measures of central tendency, i.e. the sample mean and median,
are both shown to be larger than the estimated parametric mode in Figure 1.
The sample mean, which equals 92.6 billion dollars, is obviously not a good
measure of central tendency for most large banks and savings institutions in
the United States. In particular, 40 of the 50 banks and savings institutions
in our dataset had deposits less than 92.6 billion dollars on July 2, 2010.
The sample median for this data is 40.5 billion dollars, indicating that
$50\%$ of banks in the dataset had deposits larger than 40.5 billion dollars
while the other half had deposits smaller than 40.5 billion dollars. In spite
of its high interpretability, the (sample) median is usually difficult to
visualize either from a density plot or a histogram. In contrast, it is much
easier for data analysts to locate and interpret the mode than the mean or
median in Figure 1. The estimated mode using the DTP-Student-$t$ distribution
is where the density plot reaches its peak, and is close to where the
histogram reaches to its peak. More specifically, the posterior mean of the
mode is $20$ billion dollars, suggesting that banks in the United States are
most likely to have deposits of around $20$ billion dollars during that time.

Figure 1: Deposits (in billions of dollars) of 50 banks and savings
institutions in the United States on July 2, 2010. The solid black curve is
the estimated density of the DTP-Student-$t$ distribution. The three vertical
lines mark locations of the sample mean (blue solid line), the sample median
(orange dot-dashed line), and the estimated mode (red dashed line),
respectively.

###  2.2 Modal versus mean and median regression for analyzing murder rates

As a second motivating example, we analyze a dataset from Agresti et al.,
(2021) containing the murder rate, percentage of college education, poverty
percentage and metropolitan rate for the 50 states in the United States and
the District of Columbia (D.C.) from 2003. The murder rate is defined as the
annual number of murders per $100,000$ people in the population. The poverty
percentage is the percentage of residents with income below the poverty level,
and the metropolitan rate is defined as the percentage of population living in
the metropolitan area.

At the stage of exploratory data analysis, we plotted the conditional scatter
plot matrix of the U.S. crime data in Figure 2. From the first row of the
conditional scatter plot matrix, a positive association between the poverty
percentage and the murder rate, and also a positive association between the
metropolitan rate and the murder rate are relatively evident; but an
association between the college percentage and the murder rate is harder to
perceive in the figure. Figure 2 also brings up a clear outlier, which is D.C.
Without this outlier, there appears to exist a negative association between
the college percentage and the murder rate.

Figure 2: The conditional scatter plot matrices of the U.S. crime data.

To formally investigate the association between the murder rate ($Y$) and the
aforementioned variables, we fit the following models to the U.S. crime
dataset:

$\mathbb{M}(Y\mid\boldsymbol{\beta})=\beta_{0}+\beta_{1}\times\operatorname{college}+\beta_{2}\times\operatorname{poverty}+\beta_{3}\times\operatorname{metropolitan},$

where $\mathbb{M}(\cdot)$ generically refers to the conditional mean, median,
or mode. Table 1 presents the inference results from mean/median/modal
regression models. All three models share some conclusions in common. Namely,
all models determined that there were positive associations between the
poverty percentage and the murder rate, and between the metropolitan rate and
the murder rate. However, with a posterior credible interval (CI) of
$(0.20,0.74)$, the mean regression model (specified by (5.1)) implies that
there exists a positive association between the college percentage and the
crime rate, conditionally on the other covariates in the model. We believe
that this inference result is difficult to justify, in light of existing
results from the criminology literature that conclude a negative association
between higher education attainment and crime (Lochner,, 2020; Hjalmarsson and
Lochner,, 2012). On the other hand, with a CI of $(-0.27,0.05)$, the Bayesian
median regression model (formulated in (5.2)) concludes that the college
percentage is not significantly associated with the murder rate, conditionally
on the other covariates. Our Bayesian modal regression model with the Two-
Piece scale-Student-$t$, or TPSC-Student-$t$, distribution (to be introduced
in Section 3) draws a different conclusion. With a CI of $(-0.33,-0.06)$, our
Bayesian modal regression model concludes that there is a negative association
between the college percentage and the murder rate, which is more consistent
with findings from the criminology literature. Lastly, according to the model
criterion referred to as the expected log predictive density (ELPD, introduced
in Section 4.3) in Table 1, the modal regression model based on the TPSC-
Student-$t$ likelihood yields the highest value of ELPD, indicating a better
fit for the data than the mean and median regression models.

We repeated the above analyses after removing the D.C. outlier from the data.
Now the median and modal regression models do in fact suggest a negative
association between the college percentage and the murder rate, whereas the
mean regression model insists on lack of significant association between them.
This exercise demonstrates that modal regression based on the proposed GUD
family can be even more robust to outliers than median regression, and has a
stronger potential in drawing reliable inferences and unveiling important
features of data even in the presence of extreme outliers.

Table 1:  Estimates of covariate effects for the mean/median/modal regression models fit to the U.S. crime dataset. The mean, $5\%$ quantile, and $95\%$ quantile of the posterior distribution of each covariate effect are listed under Mean, q5, and q95, respectively. ELPD stands for expected log predictive density. Regression model | ELPD | Parameter (covariate) | Mean | q5 | q95
---|---|---|---|---|---
|  |  $\beta_{1}$ (college) | 0.47 | 0.20 | 0.74
|  |  $\beta_{2}$ (poverty) | 1.14 | 0.76 | 1.51
Mean regression | -162.59 |  $\beta_{3}$ (metropolitan) | 0.07 | 0.01 | 0.12
|  |  $\beta_{1}$ (college) | -0.12 | -0.27 | 0.05
|  |  $\beta_{2}$ (poverty) | 0.44 | 0.22 | 0.67
Median regression | -133.12 |  $\beta_{3}$ (metropolitan) | 0.06 | 0.03 | 0.08
|  |  $\beta_{1}$ (college) | -0.20 | -0.33 | -0.06
|  |  $\beta_{2}$ (poverty) | 0.24 | 0.01 | 0.46
Modal regression | -123.27 |  $\beta_{3}$ (metropolitan) | 0.06 | 0.04 | 0.09

##  3 The family of general unimodal distributions

Having motivated our Bayesian modal regression framework and demonstrated its
benefits on two real-life applications in Section 2, we now formally introduce
the GUD family for Bayesian modal regression.

The probability density function (pdf) of a member of the GUD family is in the
form a mixture of two pdfs, $f_{1}$ and $f_{2}$, given by

$f(y\mid
w,\theta,\boldsymbol{\sigma}_{1},\boldsymbol{\sigma}_{2})=wf_{1}(y\mid\theta,\boldsymbol{\sigma}_{1})+(1-w)f_{2}(y\mid\theta,\boldsymbol{\sigma}_{2}).$
(3.1)

In the mixture pdf (3.1), $w\in[0,1]$ is the weight parameter,
$\theta\in\left(-\infty,+\infty\right)$ is the mode as the only location
parameter in (3.1), $\boldsymbol{\sigma}_{1}$ consists of parameters other
than the location parameter in the pdf
$f_{1}(\cdot\mid\theta,\boldsymbol{\sigma}_{1})$, and
$\boldsymbol{\sigma}_{2}$ is defined similarly for
$f_{2}(\cdot\mid\theta,\boldsymbol{\sigma}_{2})$. Clearly, the GUD family
belongs to the more general two-component mixture distribution family. One
feature of GUD that makes it stand out from the bigger family of two-component
mixture distributions is that the two component distributions of GUD share the
same location parameter $\theta$ as the mode, a feature that makes GUD
especially suitable for modal regression. In contrast, a two-component normal
mixture for instance, as a widely referenced member in the bigger family, can
be multimodal, and it is non-trivial to impose constraints on two normal
components to guarantee unimodality (Sitek,, 2016). Even after formulating a
unimodal normal mixture, its mode may not have an analytical form
(Behboodian,, 1970). Many other members in the more general two-component
mixture distribution family have the same pitfalls.

Besides unimodality, we reiterate and complement the following three
restrictions on (3.1) to make the GUD family suitable and convenient for modal
regression:

  1. (R1)

The pdfs $f_{1}(\cdot\mid\theta,\boldsymbol{\sigma}_{1})$ and
$f_{2}(\cdot\mid\theta,\boldsymbol{\sigma}_{2})$ are unimodal at $\theta$.

  2. (R2)

The pdfs $f_{1}(\cdot\mid\theta,\boldsymbol{\sigma}_{1})$ and
$f_{2}(\cdot\mid\theta,\boldsymbol{\sigma}_{2})$ are left-skewed and right-
skewed respectively.

  3. (R3)

The mixture pdf $f(\cdot\mid
w,\theta,\boldsymbol{\sigma}_{1},\boldsymbol{\sigma}_{2})$ in (3.1) is
continuous in its domain.

Restriction (R1) is already implied earlier when we stress that the two
components in (3.1) share the same location parameter $\theta$ as the finite
mode. In the context of modal regression, (R1) ensures that one can easily
link a linear predictor $\boldsymbol{X}^{\top}\boldsymbol{\beta}$ with the
conditional mode of $Y$. Because modal regression adds more value to
mean/median regression when data are skewed and contain outliers, we impose
(R2) to make members in GUD exhibit a wide range of skewness and tail
behaviors. This second restriction also solves the notorious label switching
problem that many other two-component mixture distributions suffer from,
because $f_{1}(\cdot\mid\theta,\boldsymbol{\sigma}_{1})$ and
$f_{2}(\cdot\mid\theta,\boldsymbol{\sigma}_{2})$ satisfying (R2) must come
from different distribution families in some strict sense, as opposed to, say,
both coming from the normal family. According to Theorem 1 of Teicher, (1963),
this guarantees identifiability of all parameters associated with GUD. Lastly,
(R3) eliminates ill-constructed pdfs whose mode may occur at a jump
discontinuity.

Henceforth, when a random variable $Y$ follows a distribution in the GUD
family, we state that $Y\mid
w,\theta,\boldsymbol{\sigma}_{1},\boldsymbol{\sigma}_{2}\sim\operatorname{GUD}\left(w,\theta,\boldsymbol{\sigma}_{1},\boldsymbol{\sigma}_{2}\right)$.
Like for other two-component mixture distributions, one may view
$Y=ZX_{1}+(1-Z)X_{2}$, where $X_{1}\mid\theta,\boldsymbol{\sigma}_{1}\sim
f_{1}(\cdot\mid\theta,\boldsymbol{\sigma}_{1})$,
$X_{2}\mid\theta,\boldsymbol{\sigma}_{2}\sim
f_{2}(\cdot\mid\theta,\boldsymbol{\sigma}_{2})$, and $Z\mid
w\sim\operatorname{Bernoulli}\left(w\right)$, with $Z$, $X_{1}$, and $X_{2}$
independent. This viewpoint gives rise to a data augmentation method outlined
below for generating data from a GUD effortlessly:

  1. (i)

Sample $X_{1}\mid\theta,\boldsymbol{\sigma}_{1}\sim
f_{1}\left(\cdot\mid\theta,\boldsymbol{\sigma}_{1}\right)$.

  2. (ii)

Sample $X_{2}\mid\theta,\boldsymbol{\sigma}_{2}\sim
f_{2}\left(\cdot\mid\theta,\boldsymbol{\sigma}_{2}\right)$.

  3. (iii)

Sample $Z\mid w\sim\operatorname{Bernoulli}(w)$.

  4. (iv)

$Y\leftarrow ZX_{1}+(1-Z)X_{2}$.

Having an efficient random number generation method is especially beneficial
in constructing Bayesian prediction intervals, since the most common way to
approximate the posterior predictive density is by drawing samples from the
posterior predictive distribution during the MCMC iterations. We will continue
our discussion about the Bayesian prediction intervals in Section 4.3.

Relating to existing literature, the GUD family subsumes several previously
proposed distributions, such as those introduced in Fernández and Steel,
(1998) and Rubio and Steel, (2015), as special cases. In what follows, we
detail several examples of distributions from the GUD family.

###  3.1 The flexible Gumbel distribution

For predicting extreme events, the Gumbel distribution is a popular choice in
many fields such as hydrology, earthquake forecasting, and insurance (Smith,,
2003; Vidal,, 2014; Shin et al.,, 2015). The pdf of a Gumbel distribution for
the maximum is

$f_{\operatorname{Gumbel}}\left(y\mid\theta,\sigma\right)=\frac{1}{\sigma}\exp\left\\{-\frac{y-\theta}{\sigma}-\exp\left(-\frac{y-\theta}{\sigma}\right)\right\\}\mathbb{I}\left(-\infty<y<\infty\right),$

where $\theta\in\mathbb{R}$ is the mode as the location parameter, $\sigma>0$
is the scale parameter, and $\mathbb{I}(\cdot)$ is the indicator function. To
describe data that contains a mix of extremely large and extremely small
events, Liu et al., (2022) proposed the flexible Gumbel (FG) distribution
specified by the pdf

$f_{\operatorname{FG}}\left(y\mid
w,\theta,\sigma_{1},\sigma_{2}\right)=wf_{\operatorname{Gumbel}}\left(-y\mid-\theta,\sigma_{1}\right)+(1-w)f_{\operatorname{Gumbel}}\left(y\mid\theta,\sigma_{2}\right).$
(3.2)

Mapping to (3.1), we have
$f_{1}(y|\theta,\sigma_{1})=f_{\operatorname{Gumbel}}\left(-y\mid-\theta,\sigma_{1}\right)$
as the pdf of the left-skewed Gumbel distribution for the minimum, and
$f_{2}(y|\theta,\sigma_{2})=f_{\operatorname{Gumbel}}\left(y\mid\theta,\sigma_{2}\right)$
as the pdf of the right-skewed Gumbel distribution for the maximum. We
illustrate Bayesian modal regression based on the FG likelihood in Section
6.2. The FG distribution serves as a good choice of likelihood if the data is
a mixture of extreme events, such as monthly maximum/minimum water elevation
changes, and weekly heaviest/lightest traffic on a highway.

###  3.2 The Double Two-Piece distribution

Rubio and Steel, (2015) defined the Double Two-Piece (DTP) distribution by
mixing two truncated distributions. For a pdf belonging to some location-scale
family of the form $(1/\sigma)f(\left(y-\theta\right)/\sigma\mid\delta)$ that
is unimodal at $\theta$, with a scale parameter $\sigma>0$ and a shape
parameter $\delta$, the pdf of the corresponding left $\theta$-truncated
distribution is

$f_{\operatorname{LT}}(y\mid\theta,\sigma,\delta)=\frac{2}{\sigma}f\left.\left(\frac{y-\theta}{\sigma}\right|\delta\right)\mathbb{I}(y<\theta),$
(3.3)

and the corresponding right $\theta$-truncated distribution is specified by
the following pdf,

$f_{\operatorname{RT}}(y\mid\theta,\sigma,\delta)=\frac{2}{\sigma}f\left.\left(\frac{y-\theta}{\sigma}\right|\delta\right)\mathbb{I}(y\geq\theta).$
(3.4)

By mixing the pdfs in (3.3)-(3.4), we obtain the DTP pdf as

$f_{\operatorname{DTP}}(y\mid\theta,\sigma_{1},\sigma_{2},\delta_{1},\delta_{2})=wf_{\operatorname{LT}}(y\mid\theta,\sigma_{1},\delta_{1})+(1-w)f_{\operatorname{RT}}(y\mid\theta,\sigma_{2},\delta_{2}),$
(3.5)

where

$w=\frac{\sigma_{1}f\left(0;\delta_{2}\right)}{\sigma_{1}f\left(0;\delta_{2}\right)+\sigma_{2}f\left(0;\delta_{1}\right)}$
(3.6)

as the weight chosen to produce a mixture distribution that satisfies (R3).
Restrictions (R1) and (R2) are trivially satisfied by the construction of the
left/right $\theta$-truncated pdfs in (3.3)–(3.4). Thus, DTP distributions
belong to the GUD family. Note, however, that our general GUD family (3.1)
does not require the two component densities to be truncated, as we
demonstrated earlier with the FG distribution (3.2).

As a concrete example, consider the location-scale family as the three-
parameter Student’s $t$ distributions, i.e., the non-standardized Student’s
$t$ distributions, with location parameter $\theta$, scale parameter
$\sigma>0$, and continuous degree of freedom $\delta>0$ (Geweke,, 1993).
Following (3.3) and (3.4), one has the corresponding left-skewed truncated
three-parameter Student’s $t$ distribution and the right-skewed truncated
three-parameter Student’s $t$ distribution, respectively. This leads to the
distribution defined according to (3.5) and (3.6) that we call the DTP-
Student-$t$ distribution. The DTP distribution family contains numerous
distributions, all of which are suitable for modal regression (see Rubio and
Steel, (2015) for more). In the sequel, we concentrate on the DTP-Student-$t$
distribution as a special member of the DTP distribution.

###  3.3 The Two-Piece scale distribution

By setting $\delta_{1}=\delta_{2}=\delta$ in (3.5), one obtains the pdf of a
subfamily of the DTP family proposed in Fernández and Steel, (1998), referred
to as the two-piece scale (TPSC) distribution family,

$\displaystyle f_{\operatorname{TPSC}}\left(y\mid
w,\theta,\sigma,\delta\right)$
$\displaystyle=wf_{\operatorname{LT}}\left(y\mid\theta,\sigma\sqrt{\frac{w}{1-w}},\delta\right)+$
(3.7)
$\displaystyle\phantom{=}(1-w)f_{\operatorname{RT}}\left(y\mid\theta,\sigma\sqrt{\frac{1-w}{w}},\delta\right).$

We shall point out that, in Fernández and Steel, (1998), a shape parameter
$\gamma=w^{0.5}(1-w)^{-0.5}$ is used instead of the weight parameter $w$ when
formulating the mixture pdf. We adopt the parameterization in (3.7) because we
find it more straightforward to elicit a noninformative prior for $w$ than
placing a noninformative prior on $\gamma$.

Similar to the construction of the DTP-Student-$t$ distribution, we can
construct the TPSC-Student-$t$ distribution by choosing the two component
distributions to be the left and right $\theta$-truncated three-parameter
Student’s $t$ distributions. When $w=0.5$, the TPSC-Student-$t$ distribution
converges to a normal distribution with mean $\theta$ and standard deviation
$\sigma$ as $\delta\to\infty$; and it reduces to a Cauchy distribution with
mode $\theta$ and scale parameter $\sigma$ when $\delta=1$. Hence, even as a
special case of the DTP-Student-$t$ distribution, the TPSC-Student-$t$
distribution is flexible enough to describe normally distributed data and non-
normal data with extreme outlier(s) from distributions that do not have any
finite moments. Since the TPSC-Student-$t$ has fewer parameters than the DTP-
Student-$t$ distribution, it is an adequate choice for small datasets. On the
other hand, the DTP-Student-$t$ distribution may be preferred when there is
moderate sample size. Certainly, one can conduct several rounds of modal
regression analysis assuming different unimodal distributions for the
response, such as the FG, DTP-Student-$t$, and TPSC-Student-$t$ distributions,
and then select the most appropriate model using the model selection criteria
that we introduce in Section 4.3. All of these models can be easily
implemented using the code developed for this work.

###  3.4 Type I GUD and type II GUD subfamilies

As illustrated in the preceding three subsections, the GUD family is a
generalization of several previously proposed unimodal two-component mixture
distributions. Figure 3 presents pdfs of FG, DTP-Student-$t$, and TPSC-
Student-$t$ distributions with different parameter specifications, which
encompass asymmetric, symmetric, fat-tailed, and thin-tailed densities. In
particular, the first panel in Figure 3 presents the density plot of FG
distribution with varying scale parameters of the right-skewed component. As
$\sigma_{2}$ becomes larger, the tails of FG distribution (especially the
right tail) become fatter. In the second panel, we show that the FG
distribution is symmetric if $w=0.5$ and $\sigma_{1}=\sigma_{2}$. With the
weight parameter $w$ surpassing 0.5 further, the pdf of FG distribution puts
more weight on the left-skewed part, therefore, becomes more left-skewed. In
the third panel, as the the scale parameter of the left-skewed component
$\sigma_{1}$ increases, the left tail of the DTP-Student-$t$ distribution
becomes fatter while the right tail changes little, leading to distributions
that are more left-skewed. The fourth panel shows the drastic change in the
shape of the TPSC-Student-$t$ pdf as one varies the scale parameter $\sigma$
shared by both mixture components. The last panel presents the subtle changes
in the tail behavior of TPSC-Student-$t$ distributions with different values
for the degree of freedom $\delta$ that is shared by both mixture components.

Figure 3: Density plots of different distributions in the GUD family with
different parameter specifications.

The GUD family can be further categorized into two subfamilies. Let
$\mathcal{D}_{1}$ and $\mathcal{D}_{2}$ denote the domains of
$f_{1}(\cdot\mid\theta,\boldsymbol{\sigma}_{1})$ and
$f_{2}(\cdot\mid\theta,\boldsymbol{\sigma}_{2})$ in the GUD pdf (3.1)
respectively. If $\mathcal{D}_{1}\cap\mathcal{D}_{2}\neq\varnothing$, we call
the mixture distribution the type I GUD. The FG distribution is an example of
type I GUD. In Appendix A, we present the construction of the lognormal
mixture distribution (logNM), which is another example of type I GUD. On the
other hand, if $\mathcal{D}_{1}\cap\mathcal{D}_{2}=\varnothing$, then we have
the type II GUD. The DTP distributions and the asymmetric Laplace distribution
(ALD, Koenker and Machado,, 1999) belong to this subfamily of type II GUD.
Figure 4 presents the partition of the GUD family into type I and type II GUD
subfamilies.

Figure 4: Venn diagram of the unimodal two-components mixture distributions.

##  4 Bayesian modal regression

Having defined the GUD family in Section 3, we are now in a position to
introduce our Bayesian modal regression framework. In the remainder of the
manuscript, we assume that we observe $n$ independent pairs of observations
$(\boldsymbol{X}_{1},Y_{1}),(\boldsymbol{X}_{2},Y_{2}),\ldots,(\boldsymbol{X}_{n},Y_{n})$.
Here, $\boldsymbol{X}_{i}=(X_{i1},\ldots,X_{ip})^{\top}$ denotes a vector of
$p$ covariates for the $i$th observation. We assume exchangeability in the
sense that, given
$\boldsymbol{X}=(\boldsymbol{X}_{1},\ldots,\boldsymbol{X}_{n})$ and all
parameters, $n$ observations in $\boldsymbol{Y}=(Y_{1},\ldots,Y_{n})$ are
independent. Our goal is to conduct inference about the conditional mode of
the response variable $Y$ given the covariates $\boldsymbol{X}$.

###  4.1 Prior elicitation

For all modal linear regression models in this paper, we assume that

$Y_{i}\mid\boldsymbol{X}_{i},w,\boldsymbol{\beta},\boldsymbol{\sigma}_{1},\boldsymbol{\sigma}_{2}\stackrel{{\scriptstyle\text{ind}}}{{\sim}}\operatorname{GUD}\left(w,\,\boldsymbol{X}_{i}^{\top}\boldsymbol{\beta},\,\boldsymbol{\sigma}_{1},\,\boldsymbol{\sigma}_{2}\right),\text{
for $i=1,\dots,n$,}$ (4.1)

where GUD generically refers to a member of the GUD family, and “ind” is the
acronym for “independent.” Recall that any member of the GUD family contains
the location parameter as its mode, which is
$\boldsymbol{X}_{i}^{\top}\boldsymbol{\beta}$ as the conditional mode for
$Y_{i}$ given $\boldsymbol{X}_{i}$ in (4.1).

To conduct inference for our model in (4.1), we adopt a Bayesian approach
where appropriate priors are placed on the model parameters
$(w,\boldsymbol{\beta},\boldsymbol{\sigma}_{1},\boldsymbol{\sigma}_{2})$.We
endow the weight parameter $w$ with a noninformative $\textrm{Uniform}(0,1)$
prior, and use weakly informative $\textrm{InverseGamma}(1,1)$ or
$\textrm{InverseGamma}(5,5)$ priors for all positive parameters in
$\boldsymbol{\sigma}_{1}$ and $\boldsymbol{\sigma}_{2}$. As pointed out by
Diebolt and Robert, (1994), improper priors usually lead to improper posterior
distributions for mixture distributions because of identifiability problems.
Therefore, if
$\boldsymbol{\sigma}_{1}\cap\boldsymbol{\sigma}_{2}=\varnothing$, then
improper priors should not be used for $\boldsymbol{\sigma}_{1}$ or
$\boldsymbol{\sigma}_{2}$.

On the other hand, a flat prior $p(\boldsymbol{\beta})\propto 1$ on the
regression coefficients $\boldsymbol{\beta}$ usually leads to a proper
posterior distribution because both right and left skewed components share the
same location parameter. In Section 4.2, we provide sufficient conditions
under which a flat prior can be used for $\boldsymbol{\beta}$ such that the
posterior distribution is proper. These sufficient conditions are very mild
and are likely to be satisfied in practice. All models going forward thus use
a noninformative flat prior, $p(\boldsymbol{\beta})\propto 1$, for
$\boldsymbol{\beta}$.

Revisiting the three members of the GUD family discussed in Section 3, we have
the Bayesian modal linear regression model based on the FG likelihood (3.2)
formulated as follows,

$\displaystyle
Y_{i}\mid\boldsymbol{X}_{i},w,\boldsymbol{\beta},\sigma_{1},\sigma_{2}$
$\displaystyle\stackrel{{\scriptstyle\text{ind}}}{{\sim}}\operatorname{FG}\left(w,\,\boldsymbol{X}_{i}^{\top}\boldsymbol{\beta},\,\sigma_{1},\,\sigma_{2}\right),\text{
for $i=1,\dots,n$,}$ (4.2) $\displaystyle w$
$\displaystyle\sim\operatorname{Uniform}(0,1),$
$\displaystyle\sigma_{1},\sigma_{2}$
$\displaystyle\stackrel{{\scriptstyle\text{i.i.d}}}{{\sim}}\text{InverseGamma}(1,1),$
$\displaystyle p(\boldsymbol{\beta})$ $\displaystyle\propto 1,$

where “i.i.d” refers to “independent and identically distributed.” Meanwhile,
the Bayesian modal linear regression associated with the DTP-Student-$t$
likelihood (3.5) is specified by

$\displaystyle
Y_{i}\mid\boldsymbol{X}_{i},\boldsymbol{\beta},\sigma_{1},\sigma_{2},\delta_{1},\delta_{2}$
$\displaystyle\stackrel{{\scriptstyle\text{ind}}}{{\sim}}\operatorname{DTP-
Student-}t\left(\boldsymbol{X}_{i}^{\top}\boldsymbol{\beta},\,\sigma_{1},\,\sigma_{2},\,\delta_{1},\,\delta_{2}\right),\text{
for $i=1,\dots,n$,}$ (4.3)
$\displaystyle\sigma_{1},\sigma_{2},\delta_{1},\delta_{2}$
$\displaystyle\stackrel{{\scriptstyle\text{i.i.d}}}{{\sim}}\text{InverseGamma}(1,1),$
$\displaystyle p(\boldsymbol{\beta})$ $\displaystyle\propto 1.$

Recall that the weight parameter $w$ of a DTP distribution is fully defined by
its scale and shape parameters, so in this case, there is no need to choose a
prior for $w$. Finally, the Bayesian modal linear regression associated with
the TPSC-Student-$t$ likelihood (3.7) is defined as

$\displaystyle Y_{i}\mid\boldsymbol{X}_{i},w,\boldsymbol{\beta},\sigma,\delta$
$\displaystyle\stackrel{{\scriptstyle\text{ind}}}{{\sim}}\operatorname{TPSC-
Student-}t\left(w,\,\boldsymbol{X}_{i}^{\top}\boldsymbol{\beta},\,\sigma,\,\delta\right),\text{
for $i=1,\dots,n$,}$ (4.4) $\displaystyle w$
$\displaystyle\sim\operatorname{Uniform}(0,1),$ $\displaystyle\sigma,\delta$
$\displaystyle\stackrel{{\scriptstyle\text{i.i.d}}}{{\sim}}\text{InverseGamma}(1,1),$
$\displaystyle p(\boldsymbol{\beta})$ $\displaystyle\propto 1.$

According to Proposition 4.1 in next subsection, all of the proposed Bayesian
modal regression models (4.2)–(4.4) above have proper posterior distributions.
Practitioners can construct various other Bayesian modal regression models
using the same strategy shown above. In this paper, we concentrate on the
modal regression models based on the FG, DTP-Student-$t$, and TPSC-Student-$t$
likelihoods for the sake of concreteness.

###  4.2 Sufficient conditions for posterior propriety

Since we use an improper prior, $p(\boldsymbol{\beta})\propto 1$, for the
regression coefficients $\boldsymbol{\beta}$ in our Bayesian modal regression
models, it is important to check that the posterior distribution is proper.
Theorem 4.1 gives sufficient conditions under which the GUD likelihood (3.1)
and suitably chosen priors on the model parameters lead to a proper posterior.
Proposition 4.1 extends this result to the regression setting. We stress that
our results are nonasymptotic; that is, for any fixed sample size $n$ and
fixed number of covariates $p$, the posterior is proper under very mild
conditions.

######  Theorem 4.1.

Assume that the following conditions are satisfied:

  1. (A)

Both ranges of $f_{1}\left(\cdot\mid\theta,\boldsymbol{\sigma}_{1}\right)$ and
$f_{2}\left(\cdot\mid\theta,\boldsymbol{\sigma}_{2}\right)$ are bounded from
above for any given $\theta$, $\boldsymbol{\sigma}_{1}$ and
$\boldsymbol{\sigma}_{2}$.

  2. (B)

Both $f_{1}\left(\cdot\mid\theta,\boldsymbol{\sigma}_{1}\right)$ and
$f_{2}\left(\cdot\mid\theta,\boldsymbol{\sigma}_{2}\right)$ have proper
posterior distributions under the flat prior $p(\theta)\propto 1$. That is,
$\int_{-\infty}^{+\infty}f_{1}\left(\cdot\mid\theta,\boldsymbol{\sigma}_{1}\right)d\theta<+\infty$
and
$\int_{-\infty}^{+\infty}f_{2}\left(\cdot\mid\theta,\boldsymbol{\sigma}_{2}\right)d\theta<+\infty$.

  3. (C)

The priors for $w,\boldsymbol{\sigma}_{1}$ and $\boldsymbol{\sigma}_{2}$ are
all proper.

Then, under a GUD likelihood (3.1) and a flat prior $p(\theta)\propto 1$ on
the location parameter $\theta$, the posterior distribution
$p\left(w,\theta,\boldsymbol{\sigma}_{1},\boldsymbol{\sigma}_{2}\mid
Y_{1},\dots,Y_{n}\right)$ is proper.

###### Proof.

By assumption (A), there exists a constant $c>0$ such that

$\prod_{i=1}^{n}f\left(y_{i}\mid
w,\theta,\boldsymbol{\sigma}_{1},\boldsymbol{\sigma}_{2}\right)\leq c\times
f\left(y_{1}\mid
w,\theta,\boldsymbol{\sigma}_{1},\boldsymbol{\sigma}_{2}\right).$

Applying Bayes’ theorem, we have

$\displaystyle
p\left(w,\theta,\boldsymbol{\sigma}_{1},\boldsymbol{\sigma}_{2}\mid
Y_{1},\dots,Y_{n}\right)\propto\left\\{\prod_{i=1}^{n}f\left(y_{i}\mid
w,\theta,\boldsymbol{\sigma}_{1},\boldsymbol{\sigma}_{2}\right)\right\\}p(w)p(\theta)p(\boldsymbol{\sigma}_{1})p(\boldsymbol{\sigma}_{2})$
$\displaystyle\qquad\qquad\leq
c\left\\{wf_{1}(y_{1}\mid\theta,\boldsymbol{\sigma}_{1})+(1-w)f_{2}(y_{1}\mid\theta,\boldsymbol{\sigma}_{2})\right\\}p(w)p(\theta)p(\boldsymbol{\sigma}_{1})p(\boldsymbol{\sigma}_{2}).$

By assumption (B),

$\int_{-\infty}^{+\infty}\left\\{wf_{1}(y_{1}\mid\theta,\boldsymbol{\sigma}_{1})+(1-w)f_{2}(y_{1}\mid\theta,\boldsymbol{\sigma}_{2})\right\\}p(\theta)d\theta<+\infty.$

By assumption (C), the priors for $w,\boldsymbol{\sigma}_{1}$ and
$\boldsymbol{\sigma}_{2}$ are all proper, and so the posterior
$p\left(w,\theta,\boldsymbol{\sigma}_{1},\boldsymbol{\sigma}_{2}\mid
Y_{1},\dots,Y_{n}\right)$ is proper. ∎

Condition (A) is very mild and indicated by restriction (R1) stated in Section
3 since, for $k=1,2$, $f_{k}(.|\theta,\boldsymbol{\sigma}_{k})\leq
f_{k}(\theta|\theta,\boldsymbol{\sigma}_{k})$, which is finite as implied by
the existence of a well-defined mode $\theta$ for each component distribution
in the mixture. Condition (B) can be checked easily. For example, if
$f_{2}\left(y\mid\theta,\boldsymbol{\sigma}_{2}\right)=f_{\text{Gumbel}}\left(y\mid\theta,\sigma_{2}\right)$
as in an FG distribution, one can show that
$\int_{-\infty}^{+\infty}f_{\text{Gumbel}}\left(y\mid\theta,\sigma_{2}\right)d\theta=\sigma_{2}<+\infty$.

We now extend the result in Theorem 4.1 to the regression setting. Proposition
4.1 enables us to use the noninformative flat prior
$p(\boldsymbol{\beta})\propto 1$ for the regression coefficients
$\boldsymbol{\beta}$ in Bayesian modal regression.

######  Proposition 4.1.

For independent observations
$\left(\boldsymbol{X}_{1},Y_{1}\right),\ldots,\left(\boldsymbol{X}_{n},Y_{n}\right)$,
suppose that
$\max_{1\leq i\leq n}|\boldsymbol{X}_{i}^{\top}\boldsymbol{\beta}|<+\infty$
and that the following assumptions hold:

  1. (D)

For any $i=1,\dots,n$, both ranges of
$f_{1}\left(\cdot\mid\boldsymbol{X}_{i}^{\top}\boldsymbol{\beta},\boldsymbol{\sigma}_{1}\right)$
and
$f_{2}\left(\cdot\mid\boldsymbol{X}_{i}^{\top}\boldsymbol{\beta},\boldsymbol{\sigma}_{2}\right)$
are bounded from above for any given
$\boldsymbol{X}_{i}^{\top}\boldsymbol{\beta}$, $\boldsymbol{\sigma}_{1}$ and
$\boldsymbol{\sigma}_{2}$.

  2. (E)

For any $i=1,\dots,n$, both
$f_{1}\left(\cdot\mid\boldsymbol{X}_{i}^{\top}\boldsymbol{\beta},\boldsymbol{\sigma}_{1}\right)$
and
$f_{2}\left(\cdot\mid\boldsymbol{X}_{i}^{\top}\boldsymbol{\beta},\boldsymbol{\sigma}_{2}\right)$
have proper posterior distributions under the flat prior
$p(\boldsymbol{\beta})\propto 1$. That is
$\int_{\mathbb{R}^{p}}f_{1}\left(\cdot\mid\boldsymbol{X}_{i}^{\top}\boldsymbol{\beta},\boldsymbol{\sigma}_{1}\right)d\boldsymbol{\beta}<+\infty$
and
$\int_{\mathbb{R}^{p}}f_{2}\left(\cdot\mid\boldsymbol{X}_{i}^{\top}\boldsymbol{\beta},\boldsymbol{\sigma}_{2}\right)d\boldsymbol{\beta}<+\infty$.

  3. (F)

The priors for $w,\boldsymbol{\sigma}_{1}$, and $\boldsymbol{\sigma}_{2}$ are
proper.

Then, under the GUD likelihood (3.1) and a flat prior
$p(\boldsymbol{\beta})\propto 1$ on $\boldsymbol{\beta}$, the posterior
distribution
$p\left(w,\boldsymbol{\beta},\boldsymbol{\sigma}_{1},\boldsymbol{\sigma}_{2}\mid\boldsymbol{X},\boldsymbol{Y}\right)$
is proper.

###### Proof.

By assumption (D), there must exist a constant $c>0$ such that

$\prod_{i=1}^{n}f\left(Y_{i}\mid
w,\boldsymbol{X}_{i}^{\top}\boldsymbol{\beta},\boldsymbol{\sigma}_{1},\boldsymbol{\sigma}_{2}\right)\leq
c\times f\left(Y_{1}\mid
w,\boldsymbol{X}_{1}^{\top}\boldsymbol{\beta},\boldsymbol{\sigma}_{1},\boldsymbol{\sigma}_{2}\right),$

where $f\left(Y_{1}\mid
w,\boldsymbol{X}_{i}^{\top}\boldsymbol{\beta},\boldsymbol{\sigma}_{1},\boldsymbol{\sigma}_{2}\right)$
is the pdf function defined in (3.1). Applying Bayes’ theorem, we have

$\displaystyle
p\left(w,\boldsymbol{\beta},\boldsymbol{\sigma}_{1},\boldsymbol{\sigma}_{2}\mid
Y_{1},\dots,Y_{n}\right)\propto\left\\{\prod_{i=1}^{n}f\left(Y_{i}\mid
w,\boldsymbol{X}_{i}^{\top}\boldsymbol{\beta},\boldsymbol{\sigma}_{1},\boldsymbol{\sigma}_{2}\right)\right\\}p(w)p(\boldsymbol{\beta})p(\boldsymbol{\sigma}_{1})p(\boldsymbol{\sigma}_{2})$
$\displaystyle\quad\quad\quad\leq
c\left\\{wf_{1}(Y_{1}\mid\boldsymbol{X}_{1}^{\top}\boldsymbol{\beta},\boldsymbol{\sigma}_{1})+(1-w)f_{2}(Y_{1}\mid\boldsymbol{X}_{1}^{\top}\boldsymbol{\beta},\boldsymbol{\sigma}_{2})\right\\}p(w)p(\boldsymbol{\beta})p(\boldsymbol{\sigma}_{1})p(\boldsymbol{\sigma}_{2}).$

Further, by assumption (E),

$\int_{\mathbb{R}^{p}}\left\\{wf_{1}(Y_{1}\mid\boldsymbol{X}_{1}^{\top}\boldsymbol{\beta},\boldsymbol{\sigma}_{1})+(1-w)f_{2}(Y_{1}\mid\boldsymbol{X}_{1}^{\top}\boldsymbol{\beta},\boldsymbol{\sigma}_{2})\right\\}p(\boldsymbol{\beta})d\boldsymbol{\beta}<+\infty.$

Since the priors for $w,\boldsymbol{\sigma}_{1}$ and $\boldsymbol{\sigma}_{2}$
are all proper by assumption (F), the posterior distribution
$p\left(w,\boldsymbol{\beta},\boldsymbol{\sigma}_{1},\boldsymbol{\sigma}_{2}\mid\boldsymbol{X},\boldsymbol{Y}\right)$
is also proper. ∎

Diebolt and Robert, (1994) have argued that improper priors should in general
not be used for Bayesian modeling of mixture distributions. We note, however,
that the reasoning of Diebolt and Robert, (1994) does not necessarily apply to
the location parameter $\theta$ (or the mode). This is because the mode
$\theta$ is shared by both left- and right-skewed components in our proposed
GUD family of distributions. Therefore, we are able to derive sufficient
conditions under which a totally noninformative flat prior $p(\theta)\propto
1$ or $p(\boldsymbol{\beta})\propto 1$ can still be used to infer the
conditional mode.

On the other hand, we recommend against using improper priors for any of the
non-location parameters (i.e. the shape/scale parameters) in Bayesian modal
regression based on the GUD family. We formalize this in Proposition 4.2
below. This proposition states that, for the GUD family, using an improper
prior for any shape/scale parameter that is not shared by both components
leads to an improper posterior distribution.

######  Proposition 4.2.

If
$\tau\in\left(\boldsymbol{\sigma}_{1}\cup\boldsymbol{\sigma}_{2}\right)\backslash\left(\boldsymbol{\sigma}_{1}\cap\boldsymbol{\sigma}_{2}\right)$,
then using an improper prior for $\tau$ will lead to an improper posterior
distribution.

###### Proof.

Recall that the pdf of the GUD family is

$f\left(y\mid
w,\theta,\boldsymbol{\sigma}_{1},\boldsymbol{\sigma}_{2}\right)=wf_{1}\left(y\mid\theta,\boldsymbol{\sigma}_{1}\right)+(1-w)f_{2}\left(y\mid\theta,\boldsymbol{\sigma}_{2}\right).$

Without loss of generality, we assume that $\tau\in\boldsymbol{\sigma}_{1}$
and $\tau\notin\boldsymbol{\sigma}_{2}$. Suppose that there is only one
observation, we have that

$\displaystyle\int_{\tau\in\Theta_{\tau}}p\left(w,\theta,\boldsymbol{\sigma}_{1},\boldsymbol{\sigma}_{2}\mid
y\right)d\tau$
$\displaystyle\geq\int_{\tau\in\Theta_{\tau}}(1-w)f_{2}\left(y\mid\theta,\boldsymbol{\sigma}_{2}\right)p\left(\tau\right)d\tau$
$\displaystyle=\infty.$

When there is more than one observation, binomial expansion of the GUD
likelihood gives $\prod_{i=1}^{n}f\left(y_{i}\mid
w,\theta,\boldsymbol{\sigma}_{1},\boldsymbol{\sigma}_{2}\right)\geq
C[f_{2}\left(y\mid\theta,\boldsymbol{\sigma}_{2}\right)]^{n}$, where $C>0$ is
free of $\boldsymbol{\sigma}_{1}$. Hence, for
$C[f_{2}\left(y\mid\theta,\boldsymbol{\sigma}_{2}\right)]^{n}p\left(\tau\right)$,
taking integration with respect to $\tau$ is still divergent. ∎

In Appendix A, we provide a specific example of Proposition 4.2 for the logNM
distribution.

###  4.3 Uncertainty quantification and model selection

Letting
$\boldsymbol{\Omega}=\left[w,\boldsymbol{\beta},\boldsymbol{\sigma}_{1},\boldsymbol{\sigma}_{2}\right]^{\top}$,
the posterior predictive distribution under our Bayesian modal regression
model is defined as

$\displaystyle
p\left(Y_{\operatorname{new}}\mid\boldsymbol{Y},\boldsymbol{X}\right)$
$\displaystyle=\int_{\Theta}p\left(Y_{\operatorname{new}}\mid\boldsymbol{\Omega},\boldsymbol{Y},\boldsymbol{X}\right)p\left(\boldsymbol{\Omega}\mid\boldsymbol{Y},\boldsymbol{X}\right)d\boldsymbol{\Omega}$
(4.5)
$\displaystyle=\int_{\Theta}p\left(Y_{\operatorname{new}}\mid\boldsymbol{\Omega},\boldsymbol{X}\right)p\left(\boldsymbol{\Omega}\mid\boldsymbol{Y},\boldsymbol{X}\right)d\boldsymbol{\Omega},$

where $\Theta$ denotes the parameter space, and the last line holds because of
the conditional independence of $Y_{\text{new}}$ and $\boldsymbol{Y}$.
Obtaining an approximation of the posterior predictive (4.5) is
computationally inexpensive. With the random number generation algorithm
outlined in Section 3 for the GUD family, one can easily draw samples from
$p(Y_{\text{new}}\mid\boldsymbol{\Omega},\boldsymbol{X})$ during each
iteration in our MCMC algorithm, and then obtain samples from the posterior
predictive distribution
$p\left(Y_{\operatorname{new}}\mid\boldsymbol{Y},\boldsymbol{X}\right)$. In
this paper, we use the hdi function in the R package HDInterval (R Core Team,,
2022; Meredith et al.,, 2018), whose inputs are random samples generated from
the posterior predictive distributions, to calculate the highest density
intervals (HDI) with a pre-specified nominal level of coverage probability. We
use 90% HDI intervals as the posterior prediction intervals for all
mean/median/modal regression models that we consider in Sections 2.2, 5, and
6.

Due to the inherent nature of the conditional mode, the HDI prediction
intervals from modal regression models will usually be narrower than those
constructed under mean or median regression models, while having the same
amount of coverage (Yao and Li,, 2014). From a statistical inference point of
view, this is a very attractive property of our Bayesian modal regression
models – we can obtain high coverage with tighter intervals. Prediction
intervals from mean or median regression can sometimes be very conservative
and contain many implausible values. We illustrate the benefits of more
efficient inference from modal regression in Sections 5 and 6.

As mentioned in Section 4, there are many different GUD likelihoods that a
practitioner can choose from in order to conduct Bayesian inference for modal
regression. We propose to use the Bayesian leave-one-out expected log
posterior density as a model selection criterion for selecting the “best” GUD
likelihood to use. The Bayesian leave-one-out expected log predictive density
is defined as

$\text{ELPD}=\sum_{i=1}^{n}\log p\left(Y_{i}\mid Y_{-i}\right),$ (4.6)

where $Y_{-i}$ represents all observations except the $i$-th observation. In
(4.6), “ELPD” stands for the theoretical expected log predictive density.
Intuitively, if a model fits the data well, its predicted value of $Y_{i}$
given $Y_{-i}$ should be close to the observed $Y_{i}$ and $p\left(Y_{i}\mid
Y_{-i}\right)$ should be large, for all $i=1,\dots,n$. Therefore, an adequate
model tends to yield a high ELPD.

We apply the Pareto-smoothed importance sampling method (PSIS) of Vehtari et
al., (2017) to obtain an estimate of ELPD. The PSIS estimation of ELPD has
been implemented in an R package loo, which is compatible with the Stan
programming language (Carpenter et al.,, 2017). We used the Stan programming
language interfaced with R to implement all regression analysis in this paper.
When fitting multiple competing models to the same dataset, the model with the
highest estimated ELPD is preferred. By a slight abuse of notation, we use
ELPD to refer to the estimated ELPD in all empirical study presented in this
paper.

##  5 Simulation studies

We now present a few simulation studies which show that our Bayesian modal
regression model is an excellent choice for modeling data that is heavily
skewed. Under our simulation settings, simulated data was either left-skewed
or right-skewed; and, in addition to the pronounced global conditional mode,
there was also a small local mode. We compared our Bayesian modal regression
models to Bayesian mean and median regression models. The Bayesian mean
regression used a normal likelihood, i.e.,

$\displaystyle Y_{i}\mid\boldsymbol{\beta},\sigma,\boldsymbol{X}_{i}$
$\displaystyle\stackrel{{\scriptstyle\text{ind}}}{{\sim}}\mathcal{N}(\boldsymbol{X}_{i}^{\top}\boldsymbol{\beta},\,\sigma),\text{
for $i=1,\dots,n$,}$ (5.1) $\displaystyle\sigma$
$\displaystyle\sim\text{InverseGamma}(1,1),$ $\displaystyle
p(\boldsymbol{\beta})$ $\displaystyle\propto 1.$

In line with the literature on parametric Bayesian quantile regression (Yu and
Moyeed,, 2001; Yu and Zhang,, 2005), we also implemented Bayesian median
regresson using the asymmetric Laplace distribution (ALD), with quantile
parameter $p=0.5$. That is, our Bayesian modal median regression model was

$\displaystyle Y_{i}\mid\boldsymbol{\beta},\sigma,\boldsymbol{X}_{i}$
$\displaystyle\stackrel{{\scriptstyle\text{ind}}}{{\sim}}\operatorname{ALD}(\boldsymbol{X}_{i}^{\top}\boldsymbol{\beta},\,\sigma,\,p=0.5),\text{
for $i=1,\dots,n$,}$ (5.2) $\displaystyle\sigma$
$\displaystyle\sim\text{InverseGamma}(5,5),$ $\displaystyle
p(\boldsymbol{\beta})$ $\displaystyle\propto 1.$

We stress that in our simulation studies, none of the likelihoods used for
mean, median, or modal regression was exactly the same as the data generating
mechanism. Therefore, all considered regression models are “wrong,” creating
particularly realistic yet challenging scenarios under which we could more
fairly compare the performance across these competing methods.

###  5.1 Left-skewed data

We generated $n=30$ observations from the model,

$Y_{i}=\beta_{0}+\beta_{1}X_{i}+\epsilon_{i},$

where $\beta_{0}=\beta_{1}=1$,
$\epsilon_{i}\overset{\text{i.i.d}}{\sim}0.05\mathcal{N}(-50,1^{2})+0.95\mathcal{N}(0,1^{2})$,
and $X_{1},\dots,X_{30}$ were generated independently from
$\text{Uniform}(0,1)$. We then fit the mean/median/modal regression models to
the simulated data. For modal regression, we fit the FG model (4.2), the DTP-
Student-$t$ model (4.3), and the TPSC-Student-$t$ model (4.4). Among the modal
regression models, we found that the TPSC-Student-$t$ model had the highest
ELPD. For the sake of brevity, we present only the results from the models fit
with the normal, ALD, and TPSC-Student-$t$ likelihoods.

In Figure 5, we provide the empirical coverage rate and the average width of
the posterior prediction intervals across $n=30$ observations under each of
mean/median/modal regression model. With the narrowest prediction interval for
the same amount of coverage, results from the modal regression model clearly
stand out in Figure 5. In addition, the modal regression model with the TPSC-
Student-$t$ likelihood had the largest ELPD. Therefore, it was the most
appropriate model for the simulated data among the three candidate models in
this replication.

We repeated this experiment $300$ times. Table 2 shows the mean coverage rate,
prediction interval width, and ELPD across the 300 replications. The mean
regression model with the normal likelihood had the lowest average coverage
rate and the widest posterior prediction intervals. Both the median and modal
regression models had almost identical average coverage rate. However, the
modal regression model had, on average, the narrowest prediction intervals.
Since the modal regression model with the TPSC-Student-$t$ likelihood had the
largest average ELPD, we conclude that the modal regression model based on the
TPSC-Student-$t$ provided the best model fit.

Figure 5: The gray shaded areas show the $90\%$ posterior prediction intervals for the simulated left-skewed data. The solid red line is the estimated median from the posterior predictive distribution. The prediction intervals are narrower for Bayesian modal regression.  Table 2: Comparison of Bayesian mean, median, and modal regression models fitted to left-skewed data. Results were averaged across 300 Monte-Carlo replicates of left-skewed datasets. The empirical standard error associated with each Monte-Carlo average is provided in parenthesis following the average. Likelihood (regression model) | Coverage Rate (%) | Width | ELPD
---|---|---|---
Normal (mean regression) | 93.50 (0.19) | 32.25 (1.08) | -104.76 (2.00)
ALD (median regression) | 94.69 (0.22) | 14.70 (0.49) | -84.81 (1.39)
TPSC-Student-$t$ (modal regression) | 94.70 (0.22) | 8.36 (0.21) |  -59.93 (0.64)

###  5.2 Right-skewed data

In Section 5.1, we demonstrated the advantages of Bayesian modal regression
models when the data was left-skewed. In this section, we investigate our
model’s ability to detect right-skewness. We followed the same simulation
settings as those in Section 5.1, except the residual error was
$\epsilon_{i}\overset{\text{i.i.d}}{\sim}0.025\mathcal{N}(-25,1^{2})+0.95\mathcal{N}(0,1^{2})+0.025\mathcal{N}(50,1^{2})$
so that the simulated data was right-skewed, potentially with extremely large
outliers and some outliers on the lower tail. We fit the mean/median/modal
regression models to this simulated data.

Figure 6 shows that all three models achieved a coverage rate of $93\%$ in one
experiment. However, the modal regression model with the TPSC-Student-$t$
likelihood (3.7) had the narrowest posterior prediction interval and the
largest ELPD. Table 3 shows our results averaged across $300$ repeated
experiments. The modal regression model with the TPSC-Student-$t$ likelihood
had the highest average coverage rate, the narrowest posterior prediction
intervals on average, and the largest average ELPD. Therefore, we conclude
that the Bayesian modal regression model had the best performance in this
right-skewed simulation study.

Figure 6: The gray shaded areas show the $90\%$ posterior prediction intervals for the simulated right-skewed data. The solid red line is the estimated median from the posterior predictive distribution. The prediction intervals are narrower for Bayesian modal regression.  Table 3: Comparison of Bayesian mean, median, and modal regression models fitted to right-skewed data. Results were averaged across 300 Monte-Carlo replicates of right-skewed datasets. The empirical standard error associated with each Monte-Carlo average is provided in parenthesis following the average. Likelihood (regression model) | Coverage Rate (%) | Width | ELPD
---|---|---|---
Normal (mean regression) | 93.67 (0.18) | 26.14 (0.95) | -99.41 (1.90)
ALD (median regression) | 94.71 (0.22) | 12.42 (0.40) | -79.25 (1.29)
TPSC-Student-$t$ (modal regression) | 94.73 (0.22) | 8.20 (0.20) |  -60.04 (0.65)

##  6 More data applications of Bayesian modal regression

###  6.1 Uncertainty quantification of air pollution

Inhalable particulate matter referred to as PM10 is any inhalable particle
with a diameter of 10 micrometers and smaller. PM10 includes smoke, dust, and
metals. In the environmental science literature, it has been found that low
(resp. high) wind speeds are associated with high (resp. low) PM10 values
(Cichowicz et al.,, 2020). In addition, many studies have found that PM10 is
associated with heart disease, respiratory disease, and premature death
(Schwartz,, 1999; Zhao et al.,, 2017). We analyzed the PM10 dataset from
http://lib.stat.cmu.edu/datasets/. This dataset consists of air pollution
information collected by the Norwegian Public Roads Administration at Alnabru
in Oslo, Norway from October 2001 and August 2003. The response variable is
the hourly measurements of the concentration of PM10, while the predictor is
the hourly wind speed in meters per second. We fit the following
mean/median/modal regression models to this dataset:

$\mathbb{M}(Y\mid\boldsymbol{\beta})=\beta_{0}+\beta_{1}\times\operatorname{windspeed}.$

Table 4 presents the parameter estimation results of three models. The mean
regression model implies that the association between the PM10 concentration
and the wind speed is not significant (CI of $(-2.64,0.12)$). On the other
hand, both the median and modal regression models capture the negative
association between the PM10 concentration and the wind speed (CIs of
$(-2.77,-1.03)$ and $(-1.53,-0.54)$, respectively).

Table 4: Parameter estimates obtained from the mean/median/modal regression models fitted to the air pollution data. The mean, $5\%$ quantile, and $95\%$ quantile of the posterior distribution of each regression coefficient are listed under Mean, q5, and q95, respectively. Likelihood (regression model) | Parameter | Mean | q5 | q95
---|---|---|---|---
| $\beta_{0}\left(\operatorname{intercept}\right)$ | 41.94 | 36.80 | 47.10
Normal (mean regression) | $\beta_{1}\left(\operatorname{windspeed}\right)$ | -1.27 | -2.64 | 0.12
| $\beta_{0}\left(\operatorname{intercept}\right)$ | 32.75 | 29.50 | 36.23
ALD (median regression) | $\beta_{1}\left(\operatorname{windspeed}\right)$ | -1.87 | -2.77 | -1.03
| $\beta_{0}\left(\operatorname{intercept}\right)$ | 9.67 | 7.47 | 11.75
TPSC-Student-$t$ (modal regression) | $\beta_{1}\left(\operatorname{windspeed}\right)$ | -1.01 | -1.53 | -0.54
Figure 7:  The shaded regions show the $90\%$ posterior prediction intervals
for the air pollution data. The blue dashed line is the reference line as the
minimum possible value of PM10 (zero). The solid red line represents the
estimated median from the posterior predictive distribution.

Although both median and modal regression detected the negative association
between wind speed and PM10 concentration, we see from Figure 7 that the
$90\%$ prediction intervals for median regression (and mean regression)
contain many negative values. As the minimum possible measure of PM10
concentration is zero, it is difficult to justify posterior prediction
intervals for PM10 that contain many negative values. On the other hand, the
$90\%$ prediction intervals from the modal regression model only contain a
tiny portion of negative values at very high wind speeds. This suggests that
uncertainty quantification under Bayesian modal regression is more reliable
and yields more practically meaningful results in this particular example.

###  6.2 Detecting a quadratic relationship in serum data

Table 5:  Parameter estimates from the mean/median/modal regression models fitted to the serum data. The mean, $5\%$ quantile, and $95\%$ quantile of the posterior distribution of each regression coefficient are listed under Mean, q5, and q95, respectively. Likelihood (regression model) | ELPD | Parameter | Mean | q5 | q95
---|---|---|---|---|---
|  | $\beta_{0}\left(\operatorname{intercept}\right)$ | 3.08 | 2.45 | 3.72
|  | $\beta_{1}\left(\operatorname{Age}\right)$ | 0.97 | 0.45 | 1.49
Normal (mean regression) | -627.13 | $\beta_{2}\left(\operatorname{Age}^{2}\right)$ | -0.05 | -0.13 | 0.04
|  | $\beta_{0}\left(\operatorname{intercept}\right)$ | 2.82 | 2.12 | 3.57
|  | $\beta_{1}\left(\operatorname{Age}\right)$ | 1.12 | 0.53 | 1.68
ALD (median regression) | -638.16 | $\beta_{2}\left(\operatorname{Age}^{2}\right)$ | -0.07 | -0.16 | 0.03
|  | $\beta_{0}\left(\operatorname{intercept}\right)$ | 2.37 | 1.84 | 2.89
|  | $\beta_{1}\left(\operatorname{Age}\right)$ | 1.16 | 0.72 | 1.59
FG (modal regression) | -623.15 | $\beta_{2}\left(\operatorname{Age}^{2}\right)$ | -0.11 | -0.18 | -0.03

Isaacs et al., (1983) analyzed the relationship between serum concentration
(grams per litre) of immunoglobulin-G (IgG) in 298 children aged from 6 months
to 6 years. IgG is an antibody that plays an important role in humoral and
protective immunity (van de Bovenkamp et al.,, 2016). There are ethical
difficulties in taking repeated blood samples from healthy subjects.
Therefore, researchers often use age as a proxy for determining the reference
ranges for IgG in childhood. Previously, Yu and Moyeed, (2001) analyzed serum
data and modeled IgG concentration with a quadratic model in age. In the
spirit of Yu and Moyeed, (2001), we fit the following mean/median/modal
regression models to this dataset:

$\mathbb{M}(Y\mid\boldsymbol{\beta})=\beta_{0}+\beta_{1}\times\operatorname{Age}+\beta_{2}\times\operatorname{Age}^{2}.$

Table 5 shows the parameter estimates from the models that we fit to this
data. Based on the CIs for $\beta_{2}$, we see that only the modal regression
model is able to detect the quadratic term (CI of $(-0.18,-0.03)$ for modal
regression). This finding is somewhat consistent with Royston and Altman,
(1994) who concluded that a simple linear regression model was inadequate for
this same dataset. Isaacs et al., (1983) also suggested that there was a
quadratic relationship between the square root of IgG concentration and
children’s age.

Table 5 shows that the ELPD of the modal regression model based on the FG
likelihood (4.2) is larger than the ELPD for both the mean or median
regression models. In this example, the modal regression model not only
provides a different viewpoint (i.e. that there exists a significant quadratic
relationship between IgG and age), but it also fits the dataset better
according to our model selection criterion.

##  7 Discussion

In this paper, we have introduced a unifying Bayesian modal regression
framework. Namely, we proposed a simple and flexible unimodal distribution
family called the GUD family that is suitable for Bayesian modal regression.
All members of the GUD family have a location parameter that is also the mode.
Members of this family can be either symmetric or asymmetric, either thin-
tailed or heavy-tailed, depending on values of the shape and scale parameters.

Compared to mean and quantile regression, work on Bayesian modal regression
analysis is quite scarce. Our paper aims to promote Bayesian modal regression
as a complement to these other analyses. We demonstrated that our modeling
framework based on the GUD family is very versatile and has wide applications
in many fields such as economics (the bank deposit data in Section 2.1),
criminology (the murder rate data in Section 2.2), environmental science (the
air pollution data in Section 6.1), and molecular biology (the serum data in
Section 6.2). In particular, we showed that Bayesian modal regression can
reveal structures and detect potentially significant covariate effects that
are missed by other Bayesian regression models.

To conduct Bayesian inference of the conditional mode, we provided prior
elicitation procedures, along with the sufficient conditions under which a
flat prior $p(\boldsymbol{\beta})$ on the regression coefficients
$\boldsymbol{\beta}$ can be used. We proposed a method for constructing
posterior prediction intervals and a model selection criterion based on the
posterior predictive distribution. We demonstrated that our modal regression
models provide very tight prediction intervals with high coverage, are robust
to outliers, and have excellent interpretability. Our modal regression model
framework is an especially appealing choice when the data is skewed and(or)
contains (extreme) outliers.

The modal regression models in this paper contain parametric assumptions, both
about the data likelihood and the linear relationship between the covariates
and the conditional mode. Instead of using the fully parametric models
presented in this manuscript, one may prefer to use Bayesian semiparametric
modal regression models instead. A Bayesian semiparametric modal regression
model can be constructed either by modeling the conditional mode with a
Gaussian process (i.e. we can relax the linearity assumption) and/or by
replacing the GUD likelihood with a carefully constructed infinite mixture
model that is indexed by the mode (i.e. we can relax the assumption of a known
residual error distribution). These exciting extensions to Bayesian modal
regression are the topics of ongoing work.

## Appendix

##  Appendix A The Lognormal Mixture Distribution

To demonstrate how researchers can add new members to the GUD family, we
present the construction of the lognormal mixture distribution (logNM) below.
Here, we pick the lognormal distribution because the lognormal distribution is
right-skewed and unimodal. We construct a location-shift lognormal
distribution such that the transformed lognormal distribution is still right-
skewed but has a mode at $0$. Next, we flip the location-shift lognormal
distribution at $0$ to get a left-skewed unimodal lognormal distribution.
Finally, we mix the left- and right-skewed lognormal distribution together to
construct the logNM distribution. More details of the construction of logNM
can be found in (A.1)-(A.2).

The pdf of the lognormal distribution is

$f_{\operatorname{logN}}(y\mid\mu,\nu)=\frac{1}{y\nu\sqrt{2\pi}}\exp\left(-\frac{\left(\ln(y)-\mu\right)^{2}}{2\nu^{2}}\right)\mathbb{I}(y>0),$
(A.1)

where $\mu\in(-\infty,+\infty)$ and $\nu>0$ are two parametrers, and the mode
is given by $\exp\left(\mu-\nu^{2}\right)$. We define the pdf of logNM as a
mixture of two lognormal pdfs formulated as follows,

$\displaystyle f_{\operatorname{logNM}}(y\mid
w,\theta,\boldsymbol{\sigma}_{1},\boldsymbol{\sigma}_{2})$
$\displaystyle=wf_{1}(y|\theta,\boldsymbol{\sigma}_{1})+\left(1-w\right)f_{2}(y|\theta,\boldsymbol{\sigma}_{2}),$
(A.2) $\displaystyle f_{1}(y|\theta,\boldsymbol{\sigma}_{1})$
$\displaystyle=f_{\operatorname{logN}}\left(\exp\left(\mu_{1}-\nu^{2}_{1}\right)-\left(y-\theta\right)\mid\mu_{1},\nu_{1}\right),$
$\displaystyle f_{2}(y|\theta,\boldsymbol{\sigma}_{2})$
$\displaystyle=f_{\operatorname{logN}}\left(\exp\left(\mu_{2}-\nu^{2}_{2}\right)+\left(y-\theta\right)\mid\mu_{2},\nu_{2}\right),$

where $\boldsymbol{\sigma}_{1}=\left[\mu_{1},\nu_{1}\right]^{\top}$ and
$\boldsymbol{\sigma}_{2}=\left[\mu_{2},\nu_{2}\right]^{\top}$. It can be shown
that both component distributions in (A.2) are unimodal at $\theta$, with
continuous densities over the real line, one left-skewed and the other right-
skewed. Moreover, the pdfs of the individual lognormal mixture components in
(A.2) have $0$ at the right and left boundaries of their supports. Hence, the
pdf of the logNM distribution is continuous in all of $\mathbb{R}$. Having
verified that all three restrictions (R1)-(R3) for the GUD family (introduced
in Section 3 of the main manuscript), we can proceed to use the logNM
likelihood (A.3) for Bayesian modal regression.

Figure 8 demonstrates that the logNM distributions can be asymmetric or
symmetric given different combinations of parameter values. The top panel
shows that, with an increase of $\mu_{2}$, the right tail of the logNM
distribution becomes heavier while its left tail remains almost the same. The
bottom panel shows how $\nu_{1}$ influences the amount of skewness of the
logNM distribution when all three logNM distributions are left-skewed.

Figure 8:  Density plots of the logNM distribution given different
combinations of parameter values.

Practitioners can build a Bayesian modal linear regression model based on the
logNM likelihood (A.3) as follows:

$\displaystyle
Y_{i}\mid\boldsymbol{X}_{i},w,\boldsymbol{\beta},\boldsymbol{\sigma}_{1},\boldsymbol{\sigma}_{2}$
$\displaystyle\stackrel{{\scriptstyle\text{ind}}}{{\sim}}\operatorname{logNM}\left(w,\boldsymbol{X}^{\top}_{i}\boldsymbol{\beta},\mu_{1},\nu_{1},\mu_{2},\nu_{2}\right),$
(A.3) $\displaystyle w$ $\displaystyle\sim\operatorname{Uniform}(0,1),$
$\displaystyle\nu_{1},\nu_{2}$
$\displaystyle\stackrel{{\scriptstyle\text{i.i.d}}}{{\sim}}\text{InverseGamma}(1,1),$
$\displaystyle\mu_{1},\mu_{2}$
$\displaystyle\stackrel{{\scriptstyle\text{i.i.d}}}{{\sim}}\mathcal{N}(0,100^{2}),$
$\displaystyle p(\boldsymbol{\beta})$
$\displaystyle\propto\mathcal{N}_{p}\left(\boldsymbol{0},10^{2}\times\boldsymbol{I}_{p\times
p}\right),$

where $\boldsymbol{I}_{p\times p}$ stands for the $p$ by $p$ identity matrix
and $\boldsymbol{\beta}$ is a $p$-dimension random vector. If one wishes to
use the flat prior $p(\boldsymbol{\beta})\propto 1$, then it is necessary to
check conditions (E), (F) and (G) of Lemma 4.2 in the main article.

If researchers have another right-skewed or a left-skewed distribution to work
with, then they can mimic the construction of the logNM above to propose a
member of the GUD family that works for their applications. For example, one
can use the reparameterized unimodal right-skewed Gamma distribution from
Bourguignon et al., (2020) to construct a new type I GUD and a corresponding
modal regression model.

In the following lemma, we show why flat priors cannot be used for $\mu_{1}$
and(or) $\mu_{2}$ in the Bayesian modal regression model based on the logNM
distribution. This simple example demonstrates why we should avoid placing
improper priors on the non-location parameters in Bayesian modal regression
models based on the GUD family.

######  Proposition A.1.

Endowing $\mu_{1}$ and(or) $\mu_{2}$ with flat priors $p(\mu_{1})\propto 1$
and(or) $p(\mu_{2})\propto 1$ leads to an improper posterior distribution
under the logNM model (A.3).

###### Proof.

We want to show that

$\int_{-\infty}^{+\infty}\prod_{i=1}^{n}f_{\operatorname{logNM}}\left(y_{i}\mid
w,\theta,\sigma_{1},\sigma_{2},\mu_{1},\mu_{2}\right)d\mu_{1}=+\infty,$

and(or),

$\int_{-\infty}^{+\infty}\prod_{i=1}^{n}f_{\operatorname{logNM}}\left(y_{i}\mid
w,\theta,\sigma_{1},\sigma_{2},\mu_{1},\mu_{2}\right)d\mu_{2}=+\infty.$

Since

$\displaystyle f_{\operatorname{logNM}}\left(y\mid
w,\theta,\sigma_{1},\sigma_{2},\mu_{1},\mu_{2}\right)$
$\displaystyle=wf_{\operatorname{logN}}\left(\exp\left(\mu_{1}-\nu^{2}_{1}\right)-\left(y-\theta\right)\mid\theta,\mu_{1},\nu_{1}\right)+$
$\displaystyle\phantom{=}\left(1-w\right)f_{\operatorname{logN}}\left(\exp\left(\mu_{2}-\nu^{2}_{2}\right)+\left(y-\theta\right)\mid\theta,\mu_{2},\nu_{2}\right),$

and any pdf must be nonnegative, it suffices to show that

$\int_{-\infty}^{+\infty}f_{\operatorname{logN}}\left(\exp\left(\mu_{1}-\nu^{2}_{1}\right)-\left(y-\theta\right)\mid\theta,\mu_{1},\nu_{1}\right)d\mu_{2}=\infty,$

and(or)

$\int_{-\infty}^{+\infty}f_{\operatorname{logN}}\left(\exp\left(\mu_{2}-\nu^{2}_{2}\right)+\left(y-\theta\right)\mid\theta,\mu_{2},\nu_{2}\right)d\mu_{1}=\infty.$

Both the integrals above are non-finite. This completes the proof. ∎

Following the same arguments as those in Proposition A.1, one can show that
improper priors such as $p(\nu_{1})\propto 1/\nu_{1}$ and(or)
$p(\nu_{2})\propto 1/\nu_{2}$ will also lead to an improper posterior
distribution. As stated in Proposition 4.2 in the main article, a general rule
is that, for the Bayesian modal regression models based on the GUD family,
using improper prior(s) for any parameter in
$\left(\boldsymbol{\sigma}_{1}\cup\boldsymbol{\sigma}_{2}\right)\backslash\left(\boldsymbol{\sigma}_{1}\cap\boldsymbol{\sigma}_{2}\right)$
leads to an improper posterior distribution. Here, $A\backslash B=A\cap B^{c}$
denotes a collection of elements in $A$ but not in $B$.

##  Appendix B A Short Note About Markov Chain Monte Carlo (MCMC)

Readers who are familiar with Bayesian modeling of mixture distributions may
wonder why we do not use the data augmentation “trick” to design a specific
MCMC algorithm for modal regression models based on the GUD family. The
problem lies in the type II distributions of the GUD family. If
$\mathcal{D}_{1}\cap\mathcal{D}_{2}=\varnothing$, then the latent variable
conditional on other parameters and observed data becomes a degenerate random
variable. This degenerate random variable will not behave randomly.

To demonstrate that we have a degenerate random variable, let us consider a
simple case with a single observation. Recall that the type II GUD has the pdf

$f\left(y\mid
w,\theta,\boldsymbol{\sigma}_{1},\boldsymbol{\sigma}_{2}\right)=wf_{1}\left(y\mid\theta,\boldsymbol{\sigma}_{1}\right)I\left(y<\theta\right)+(1-w)f_{2}\left(y\mid\theta,\boldsymbol{\sigma}_{2}\right)I\left(y\geq\theta\right).$

Introducing the latent variable $z$, we have the joint pdf as

$f\left(y,z\mid
w,\theta,\boldsymbol{\sigma}_{1},\boldsymbol{\sigma}_{2}\right)=\left[wf_{1}\left(y\mid\theta,\boldsymbol{\sigma}_{1}\right)I\left(y<\theta\right)\right]^{z}\left[(1-w)f_{2}\left(y\mid\theta,\boldsymbol{\sigma}_{2}\right)I\left(y\geq\theta\right)\right]^{1-z}.$

The conditional distribution of $z$ is then

$\displaystyle p\left(z\mid
w,\theta,\boldsymbol{\sigma}_{1},\boldsymbol{\sigma}_{2},y\right)$
$\displaystyle\propto f\left(y,z\mid
w,\theta,\boldsymbol{\sigma}_{1},\boldsymbol{\sigma}_{2}\right)$
$\displaystyle\sim\operatorname{Bernoulli}\left(\frac{wf_{1}\left(y\mid\theta,\boldsymbol{\sigma}_{1}\right)I\left(y<\theta\right)}{wf_{1}\left(y\mid\theta,\boldsymbol{\sigma}_{1}\right)I\left(y<\theta\right)+(1-w)f_{2}\left(y\mid\theta,\boldsymbol{\sigma}_{2}\right)I\left(y\geq\theta\right)}\right).$

Similarly, the conditional distribution of $\theta$ is

$p\left(\theta\mid
w,\boldsymbol{\sigma}_{1},\boldsymbol{\sigma}_{2},y,z\right)\propto
f\left(y,z\mid
w,\theta,\boldsymbol{\sigma}_{1},\boldsymbol{\sigma}_{2}\right)p\left(\theta\right).$

The conditional mean of the latent variable $z$ can only be 0 or 1 since
$y<\theta$ and $y\geq\theta$ cannot happen at the same time. Hence, the latent
variable becomes a degenerate random variable during the MCMC iterations.
Without loss of generality, let us assume $z=1$ in the first iteration of the
MCMC. It is not hard to see that, during the MCMC iterations, the updated
values of $\theta$ can only be larger than $y$. This leads to $z$ being equal
to 1 for the rest of the iterations in the MCMC algorithm. However, if the
true value of $\theta$ is smaller than $y$, then the MCMC chain will never
reach the true value. Similarly, if $z=0$ in the first iteration of the MCMC,
then all updated values of $\theta$ can only be smaller than $y$.

In conclusion, for Bayesian modal regression models based on the type II GUD
subfamily, the latent variable data augmentation algorithm will not explore
the whole parameter space. Hence, it is non-ergodic. As a consequence of this,
we do not use the data augmentation “trick” for mixture models in our MCMC
algorithm. Instead, we use the No-U-Turn Sampler implemented in the Stan
software (Hoffman et al.,, 2014; Carpenter et al.,, 2017).

##  Appendix C Convergence Diagnostics for the Real Data Applications and
Simulation Studies

In this section, we include more details about posterior inference,
convergence diagnostics, and traceplots for the four data application examples
and two simulation studies from the main manuscript. The rhat, which has the
theoretical minimum value as 1, is a statistic measuring the convergence of
the MCMC chains. To obtain reliable posterior inference, it is recommended
that rhat should be near 1 or at least less than 1.1 (page 287 of Gelman et
al., (2013)). The ess_bulk and ess_tail are the bulk and tail effective sample
size respectively. The ess_tail is defined as the minimum of the effective
sample sizes for the 5% and 95% quantiles. The recommended lower threshold for
ess_bulk and ess_tail is 400 (Vehtari et al.,, 2021). All of the rhat,
ess_bulk and ess_tail summary statistics from our analyses in the main
manuscript satisfy these recommended thresholds.

All of the traceplots presented in this section also indicate that the MCMC
chains mixed well, with no convergence issues. We ran four MCMC chains for
each model that was fit and used the combined MCMC samples to approximate the
posterior distributions. For the intercept-only regression model fit to the
bank deposits data (Section 2.1 of the main manuscript), we set the number of
warmup iterations as $10{,}000$ and the number of post-warmup iterations as
$20{,}000$ for each chain. For all other models in the main article, we set
the number of warmup iterations as $10{,}000$ and the number of post-warmup
iterations as $10{,}000$ for each of the four MCMC chains.

###  C.1 Bank Deposits Application from Section 2.1

Figure 9:  Traceplots for the modal regression model based on the DTP-
Student-$t$ likelihood fit to the bank deposits data.

    Ψ# A tibble: 6  10
    Ψvariable     mean   median    sd   mad       q5     q95  rhat ess_bulk ess_tail
    Ψ<chr>       <dbl>    <dbl> <dbl> <dbl>    <dbl>   <dbl> <dbl>    <dbl>    <dbl>
    Ψ1 lp__     -288.    -288.    1.63  1.43  -291.    -286.    1.00   36112.   47263.
    Ψ2 sigma1      1.20     1.08  0.559 0.437    0.567    2.25  1.00   51221.   45698.
    Ψ3 sigma2     20.5     20.0   5.37  5.20    12.6     30.1   1.00   46745.   49890.
    Ψ4 delta1      1.34     1.16  0.745 0.503    0.597    2.67  1.00   70792.   47613.
    Ψ5 delta2      0.915    0.888 0.211 0.198    0.621    1.30  1.00   51342.   54038.
    Ψ6 theta      18.7     18.7   1.80  1.67    15.9     21.6   1.00   42904.   40976.

###  C.2 Crime Rate Application from Section 2.2

(a) Traceplots for the mean regression model

(b) Traceplots for the median regression model

(c) Traceplots for the modal regression model based on the TPSC-Student-$t$
likelihood

Figure 10:  Traceplots for the mean/median/modal regression models fit to the
crime data.

    Ψ## mean regression - Normal likelihood
    Ψ
    Ψ# A tibble: 4  10
    Ψvariable     mean   median     sd    mad       q5     q95  rhat ess_bulk ess_tail
    Ψ<chr>       <dbl>    <dbl>  <dbl>  <dbl>    <dbl>   <dbl> <dbl>    <dbl>    <dbl>
    Ψ1 alpha    -24.2    -24.2    5.33   5.28   -33.0    -15.4    1.00   15362.   20037.
    Ψ2 beta[1]    0.467    0.466  0.163  0.161    0.200    0.738  1.00   17845.   20507.
    Ψ3 beta[2]    1.14     1.14   0.227  0.225    0.765    1.51   1.00   18553.   22035.
    Ψ4 beta[3]    0.0677   0.0677 0.0341 0.0340   0.0121   0.124  1.00   23524.   24312.
    Ψ
    Ψ## median regression - ALD likelihood
    Ψ
    Ψ# A tibble: 4  10
    Ψvariable    mean  median     sd    mad      q5    q95  rhat ess_bulk ess_tail
    Ψ<chr>      <dbl>   <dbl>  <dbl>  <dbl>   <dbl>  <dbl> <dbl>    <dbl>    <dbl>
    Ψ1 alpha    -1.34   -1.18   3.39   3.28   -7.16   3.89    1.00    9895.   11638.
    Ψ2 beta[1]  -0.118  -0.123  0.0964 0.0933 -0.268  0.0471  1.00   10933.   13580.
    Ψ3 beta[2]   0.437   0.432  0.138  0.128   0.216  0.673   1.00   11652.   13166.
    Ψ4 beta[3]   0.0555  0.0554 0.0161 0.0158  0.0294 0.0818  1.00   17520.   18274.
    Ψ
    Ψ## modal regression - TPSC-Student-t
    Ψ
    Ψ# A tibble: 4  10
    Ψvariable    mean  median     sd    mad       q5     q95  rhat ess_bulk ess_tail
    Ψ<chr>      <dbl>   <dbl>  <dbl>  <dbl>    <dbl>   <dbl> <dbl>    <dbl>    <dbl>
    Ψ1 alpha     1.12    1.25   2.68   2.62   -3.56     5.29    1.00   13814.   16690.
    Ψ2 beta[1]  -0.199  -0.204  0.0825 0.0805 -0.327   -0.0552  1.00   14342.   17771.
    Ψ3 beta[2]   0.243   0.250  0.138  0.144   0.00688  0.459   1.00   10233.   15650.
    Ψ4 beta[3]   0.0636  0.0624 0.0153 0.0149  0.0403   0.0906  1.00   10145.    9351.

###  C.3 Air Pollution Application from Section 6

(a) Traceplots for the mean regression model

(b) Traceplots for the median regression model

(c) Traceplots for the modal regression model based on the TPSC-Student-$t$
likelihood

Figure 11:  Traceplots for the mean/median/modal regression models fit to the
air pollution data.

    Ψ## mean regression - Normal likelihood
    Ψ
    Ψ# A tibble: 2  10
    Ψvariable  mean median    sd   mad    q5    q95  rhat ess_bulk ess_tail
    Ψ<chr>    <dbl>  <dbl> <dbl> <dbl> <dbl>  <dbl> <dbl>    <dbl>    <dbl>
    Ψ1 alpha    41.9   41.9  3.12  3.14  36.8  47.1    1.00   15822.   19673.
    Ψ2 beta[1]  -1.27  -1.26 0.841 0.849 -2.64  0.116  1.00   15614.   18993.
    Ψ
    Ψ
    Ψ## median regression - ALD likelihood
    Ψ
    Ψ# A tibble: 2  10
    Ψvariable  mean median    sd   mad    q5   q95  rhat ess_bulk ess_tail
    Ψ<chr>    <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>    <dbl>    <dbl>
    Ψ1 alpha    32.8   32.7  2.05  2.03  29.5  36.2   1.00   13984.   16088.
    Ψ2 beta[1]  -1.87  -1.85 0.530 0.524 -2.77 -1.03  1.00   14043.   15873.
    Ψ
    Ψ## modal regression - TPSC-Student-t
    Ψ
    Ψ# A tibble: 2  10
    Ψvariable  mean median    sd   mad    q5    q95  rhat ess_bulk ess_tail
    Ψ<chr>    <dbl>  <dbl> <dbl> <dbl> <dbl>  <dbl> <dbl>    <dbl>    <dbl>
    Ψ1 alpha     9.67  9.76  1.44  1.18   7.46 11.7    1.00    4268.    1787.
    Ψ2 beta[1]  -1.01 -0.997 0.309 0.302 -1.53 -0.538  1.00    4581.    1832.

###  C.4 Serum Data Application from Section 6.2

(a) Traceplots for the mean regression model

(b) Traceplots for the median regression model

(c)  Traceplots for the modal regression model based on the TPSC-Student-$t$
likelihood

Figure 12:  Traceplots for the mean/median/modal regression models fit to the
serum data.

    Ψ## mean regression - Normal likelihood
    Ψ
    Ψ# A tibble: 3  10
    Ψvariable    mean  median     sd    mad     q5    q95  rhat ess_bulk ess_tail
    Ψ<chr>      <dbl>   <dbl>  <dbl>  <dbl>  <dbl>  <dbl> <dbl>    <dbl>    <dbl>
    Ψ1 alpha     3.08    3.08   0.384  0.380   2.45  3.72    1.00   10706.   15240.
    Ψ2 beta[1]   0.969   0.969  0.313  0.311   0.451 1.48    1.00    9820.   13370.
    Ψ3 beta[2]  -0.0458 -0.0460 0.0510 0.0508 -0.130 0.0385  1.00   10328.   14088.
    Ψ
    Ψ
    Ψ## median regression - ALD likelihood
    Ψ
    Ψ# A tibble: 3  10
    Ψvariable    mean  median     sd    mad     q5    q95  rhat ess_bulk ess_tail
    Ψ<chr>      <dbl>   <dbl>  <dbl>  <dbl>  <dbl>  <dbl> <dbl>    <dbl>    <dbl>
    Ψ1 alpha     2.82    2.81   0.441  0.449   2.12  3.57    1.00    8105.   12034.
    Ψ2 beta[1]   1.12    1.13   0.348  0.352   0.527 1.68    1.00    7692.   10288.
    Ψ3 beta[2]  -0.0656 -0.0668 0.0574 0.0576 -0.157 0.0314  1.00    8134.   10823.
    Ψ
    Ψ## modal regression - TPSC-Student-t
    Ψ
    Ψ# A tibble: 3  10
    Ψvariable   mean median     sd    mad     q5     q95  rhat ess_bulk ess_tail
    Ψ<chr>     <dbl>  <dbl>  <dbl>  <dbl>  <dbl>   <dbl> <dbl>    <dbl>    <dbl>
    Ψ1 alpha     2.37   2.37  0.320  0.316   1.84   2.89    1.00   12565.   17546.
    Ψ2 beta[1]   1.15   1.16  0.264  0.264   0.720  1.59    1.00   11055.   16023.
    Ψ3 beta[2]  -0.107 -0.107 0.0441 0.0438 -0.179 -0.0338  1.00   11692.   17288.

###  C.5 Left-Skewed Simulation Study from Section 5.1

(a) Traceplots for the mean regression model

(b) Traceplots for the median regression model

(c) Traceplots for the modal regression model based on the TPSC-Student-$t$
likelihood

Figure 13:  Traceplots for the mean/median/modal regression models from the
left-skewed simulation study.

    Ψ## mean regression - Normal likelihood
    Ψ
    Ψ# A tibble: 2  10
    Ψvariable  mean median    sd   mad    q5   q95  rhat ess_bulk ess_tail
    Ψ<chr>    <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>    <dbl>    <dbl>
    Ψ1 alpha    -2.32  -2.32  2.37  2.33 -6.21  1.59  1.00   36206.   27772.
    Ψ2 beta[1]   2.64   2.66  3.67  3.59 -3.36  8.66  1.00   34513.   27175.
    Ψ
    Ψ
    Ψ## median regression - ALD likelihood
    Ψ
    Ψ# A tibble: 2  10
    Ψvariable  mean median    sd   mad      q5   q95  rhat ess_bulk ess_tail
    Ψ<chr>    <dbl>  <dbl> <dbl> <dbl>   <dbl> <dbl> <dbl>    <dbl>    <dbl>
    Ψ1 alpha    0.653  0.648 0.468 0.445 -0.0964  1.43  1.00   31916.   25876.
    Ψ2 beta[1]  0.751  0.746 0.755 0.693 -0.483   2.00  1.00   34392.   25226.
    Ψ
    Ψ## modal regression - TPSC-Student-t
    Ψ
    Ψ# A tibble: 2  10
    Ψvariable  mean median    sd   mad     q5   q95  rhat ess_bulk ess_tail
    Ψ<chr>    <dbl>  <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl>    <dbl>    <dbl>
    Ψ1 alpha    0.847  0.819 0.386 0.364 0.265   1.52  1.00   20058.   19143.
    Ψ2 beta[1]  0.639  0.644 0.334 0.321 0.0800  1.18  1.00   29384.   25974.

###  C.6 Right-Skewed Simulation Study from Section 5.2

(a) Trraceplots for the mean regression model

(b) Traceplots for the median regression model

(c) Traceplots for the modal regression model based on the TPSC-Student-$t$
likelihood

Figure 14:  Traceplots for the mean/median/modal regression models from the
right-skewed simulation study.

    Ψ## mean regression - Normal likelihood
    Ψ
    Ψ# A tibble: 2  10
    Ψvariable  mean median    sd   mad    q5   q95  rhat ess_bulk ess_tail
    Ψ<chr>    <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>    <dbl>    <dbl>
    Ψ1 alpha     1.83   1.82  2.01  1.97 -1.46  5.12  1.00   37079.   27813.
    Ψ2 beta[1]   3.12   3.11  3.43  3.35 -2.51  8.76  1.00   34832.   26126.
    Ψ
    Ψ
    Ψ## median regression - ALD likelihood
    Ψ
    Ψ# A tibble: 2  10
    Ψvariable  mean median    sd   mad     q5   q95  rhat ess_bulk ess_tail
    Ψ<chr>    <dbl>  <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl>    <dbl>    <dbl>
    Ψ1 alpha    0.976  0.986 0.470 0.464  0.192  1.73  1.00   30065.   23935.
    Ψ2 beta[1]  0.473  0.487 0.944 0.925 -1.12   1.98  1.00   28958.   26648.
    Ψ
    Ψ## modal regression - TPSC-Student-t
    Ψ
    Ψ# A tibble: 2  10
    Ψvariable  mean median    sd   mad     q5   q95  rhat ess_bulk ess_tail
    Ψ<chr>    <dbl>  <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl>    <dbl>    <dbl>
    Ψ1 alpha    1.33   1.41  0.532 0.517  0.328  2.05  1.00   14412.   15244.
    Ψ2 beta[1]  0.257  0.211 0.589 0.612 -0.633  1.28  1.00   20263.   24407.

## References

  * Agresti et al., (2021) Agresti, A., Franklin, C., and Klingenberg, B. (2021).  Statistics: The Art and Science of Learning from Data.  Pearson Education, 5 edition.
  * Aristodemou, (2014) Aristodemou, K. (2014).  New regression methods for measures of central tendency.  PhD thesis.
  * Behboodian, (1970) Behboodian, J. (1970).  On the modes of a mixture of two normal distributions.  Technometrics, 12(1):131–139.
  * Benhabib and Bisin, (2018) Benhabib, J. and Bisin, A. (2018).  Skewed wealth distributions: Theory and empirics.  Journal of Economic Literature, 56(4):1261–91.
  * Bourguignon et al., (2020) Bourguignon, M., Leão, J., and Gallardo, D. I. (2020).  Parametric modal regression with varying precision.  Biometrical Journal, 62(1):202–220.
  * Carpenter et al., (2017) Carpenter, B., Gelman, A., Hoffman, M. D., Lee, D., Goodrich, B., Betancourt, M., Brubaker, M., Guo, J., Li, P., and Riddell, A. (2017).  Stan: A probabilistic programming language.  Journal of Statistical Software, 76(1).
  * Chen, (2018) Chen, Y.-C. (2018).  Modal regression using kernel density estimation: A review.  Wiley Interdisciplinary Reviews: Computational Statistics, 10(4):e1431.
  * Cichowicz et al., (2020) Cichowicz, R., Wielgosiński, G., and Fetter, W. (2020).  Effect of wind speed on the level of particulate matter pm10 concentration in atmospheric air during winter season in vicinity of large combustion plant.  Journal of Atmospheric Chemistry, 77(1):35–48.
  * Damien et al., (2017) Damien, P., Walker, S., et al. (2017).  Bayesian mode regression using mixtures of triangular densities.  Journal of Econometrics, 197(2):273–283.
  * Diebolt and Robert, (1994) Diebolt, J. and Robert, C. P. (1994).  Estimation of finite mixture distributions through Bayesian sampling.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 56(2):363–375.
  * Fernández and Steel, (1998) Fernández, C. and Steel, M. F. (1998).  On bayesian modeling of fat tails and skewness.  Journal of the American Statistical Association, 93(441):359–371.
  * Gelman et al., (2013) Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., and Rubin, D. B. (2013).  Bayesian Data Analysis.  CRC Press.
  * Geweke, (1993) Geweke, J. (1993).  Bayesian treatment of the independent student-t linear model.  Journal of applied econometrics, 8(S1):S19–S40.
  * Hjalmarsson and Lochner, (2012) Hjalmarsson, R. and Lochner, L. (2012).  The impact of education on crime: International evidence.  CESifo DICE Report, 10(2):49–55.
  * Hoffman et al., (2014) Hoffman, M. D., Gelman, A., et al. (2014).  The No-U-Turn sampler: Adaptively setting path lengths in Hamiltonian Monte Carlo.  Journal of Machine Learning Research, 15(1):1593–1623.
  * Isaacs et al., (1983) Isaacs, D., Altman, D., Tidmarsh, C., Valman, H., and Webster, A. (1983).  Serum immunoglobulin concentrations in preschool children measured by laser nephelometry: Reference ranges for igg, iga, igm.  Journal of Clinical Pathology, 36(10):1193–1196.
  * Koenker et al., (2017) Koenker, R., Chernozhukov, V., He, X., and Peng, L. (2017).  Handbook of Quantile Regression.  Chapman and Hall/CRC.
  * Koenker and Machado, (1999) Koenker, R. and Machado, J. A. (1999).  Goodness of fit and related inference processes for quantile regression.  Journal of the american statistical association, 94(448):1296–1310.
  * Lee, (1989) Lee, M.-j. (1989).  Mode regression.  Journal of Econometrics, 42(3):337–349.
  * Lee, (1993) Lee, M.-J. (1993).  Quadratic mode regression.  Journal of Econometrics, 57(1-3):1–19.
  * Liu et al., (2022) Liu, Q., Huang, X., and Zhou, H. (2022).  The flexible gumbel distribution: A new model for heavy-tailed data.  Manuscript submitted for publication.
  * Lochner, (2020) Lochner, L. (2020).  Education and crime.  In The Economics of Education, pages 109–117. Elsevier.
  * Menezes et al., (2021) Menezes, A. F., Mazucheli, J., and Chakraborty, S. (2021).  A collection of parametric modal regression models for bounded data.  Journal of Biopharmaceutical Statistics, 31(4):490–506.
  * Meredith et al., (2018) Meredith, M., Kruschke, J., and Meredith, M. M. (2018).  Package ‘hdinterval’.  Highest (Posterior) Density Intervals.
  * R Core Team, (2022) R Core Team (2022).  R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing, Vienna, Austria.
  * Royston and Altman, (1994) Royston, P. and Altman, D. G. (1994).  Regression using fractional polynomials of continuous covariates: parsimonious parametric modelling.  Journal of the Royal Statistical Society: Series C (Applied Statistics), 43(3):429–453.
  * Rubio and Steel, (2015) Rubio, F. and Steel, M. (2015).  Bayesian modelling of skewness and kurtosis with two-piece scale and shape distributions.  Electronic Journal of Statistics, 9(2):1884–1912.
  * Sager and Thisted, (1982) Sager, T. W. and Thisted, R. A. (1982).  Maximum likelihood estimation of isotonic modal regression.  The Annals of Statistics, pages 690–707.
  * Schwartz, (1999) Schwartz, J. (1999).  Air pollution and hospital admissions for heart disease in eight us counties.  Epidemiology, pages 17–22.
  * Shin et al., (2015) Shin, J. Y., Chen, S., and Kim, T.-W. (2015).  Application of bayesian markov chain monte carlo method with mixed gumbel distribution to estimate extreme magnitude of tsunamigenic earthquake.  KSCE Journal of Civil Engineering, 19(2):366–375.
  * Siegel, (2016) Siegel, A. F. (2016).  Practical Business Statistics.  Academic Press.
  * Sitek, (2016) Sitek, G. (2016).  The modes of a mixture of two normal distributions.  Silesian Journal of Pure and Applied Mathematics, 6(1):59–67.
  * Smith, (2003) Smith, R. L. (2003).  Statistics of extremes, with applications in environment, insurance, and finance.  Extreme Values in Finance, Telecommunications, and the Environment, pages 20–97.
  * Teicher, (1963) Teicher, H. (1963).  Identifiability of finite mixtures.  The Annals of Mathematical Statistics, pages 1265–1269.
  * van de Bovenkamp et al., (2016) van de Bovenkamp, F. S., Hafkenscheid, L., Rispens, T., and Rombouts, Y. (2016).  The emerging importance of igg fab glycosylation in immunity.  The Journal of Immunology, 196(4):1435–1441.
  * Vehtari et al., (2017) Vehtari, A., Gelman, A., and Gabry, J. (2017).  Practical bayesian model evaluation using leave-one-out cross-validation and waic.  Statistics and Computing, 27(5):1413–1432.
  * Vehtari et al., (2021) Vehtari, A., Gelman, A., Simpson, D., Carpenter, B., and Bürkner, P.-C. (2021).  Rank-normalization, folding, and localization: an improved r for assessing convergence of MCMC (with discussion).  Bayesian Analysis, 16(2):667–718.
  * Vidal, (2014) Vidal, I. (2014).  A bayesian analysis of the gumbel distribution: An application to extreme rainfall data.  Stochastic Environmental Research and Risk Assessment, 28(3):571–582.
  * Xiang and Yao, (2022) Xiang, S. and Yao, W. (2022).  Modal regression for skewed, truncated, or contaminated data with outliers.  In Advances and Innovations in Statistics and Data Science, pages 257–273. Springer.
  * Yao and Li, (2014) Yao, W. and Li, L. (2014).  A new regression model: Modal linear regression.  Scandinavian Journal of Statistics, 41(3):656–671.
  * Yu and Aristodemou, (2012) Yu, K. and Aristodemou, K. (2012).  Bayesian mode regression.  arXiv preprint arXiv:1208.0579.
  * Yu and Moyeed, (2001) Yu, K. and Moyeed, R. A. (2001).  Bayesian quantile regression.  Statistics & Probability Letters, 54(4):437–447.
  * Yu and Zhang, (2005) Yu, K. and Zhang, J. (2005).  A three-parameter asymmetric laplace distribution and its extension.  Communications in Statistics—Theory and Methods, 34(9-10):1867–1879.
  * Zhao et al., (2017) Zhao, Y., Cheng, Z., Lu, Y., Chang, X., Chan, C., Bai, Y., Zhang, Y., and Cheng, N. (2017).  Pm10 and pm2. 5 particles as main air pollutants contributing to rising risks of coronary heart disease: A systematic review.  Environmental Technology Reviews, 6(1):174–185.
  * Zhou and Huang, (2020) Zhou, H. and Huang, X. (2020).  Parametric mode regression for bounded responses.  Biometrical Journal, 62(7):1791–1809.
  * Zhou and Huang, (2022) Zhou, H. and Huang, X. (2022).  Bayesian beta regression for bounded responses with unknown supports.  Computational Statistics & Data Analysis, 167:107345.