# MM Algorithms for Statistical Estimation in Quantile Regression

Yifan CHENG<EMAIL_ADDRESS>Department of Statistics and Data Science,
National University of Singapore Anthony KUK Yung Cheung  Department of
Statistics and Data Science, National University of Singapore

###### Abstract

Quantile regression [28] is a robust and practically useful way to efficiently
model quantile varying correlation and predict varied response quantiles of
interest. This article constructs and tests MM algorithms, which are simple to
code and have been suggested superior to some other prominent quantile
regression methods in nonregularized problems [35], in an array of quantile
regression settings including linear (modeling different quantile coefficients
both separately and simultaneously), nonparametric, regularized, and monotone
quantile regression. Applications to various real data sets and two simulation
studies comparing MM to existing tested methods have corroborated our
algorithms’ effectiveness. We have made one key advance by generalizing our MM
algorithm to efficiently fit easy-to-predict-and-interpret parametric quantile
regression models for data sets exhibiting manifest complicated nonlinear
correlation patterns, which has not yet been covered by current literature to
the best of our knowledge.

###### Contents

  1. 1 Introduction
  2. 2 MM Algorithms
    1. 2.1 An Overview of the MM Algorithm
    2. 2.2 An MM Algorithm in Quantile Regression
  3. 3 Linear Quantile Regression
    1. 3.1 Estimating Different Quantile Coefficients Separately
    2. 3.2 Estimating Different Quantile Coefficients Simultaneously
    3. 3.3 Simulation Studies
  4. 4 Nonparametric Quantile Regression
    1. 4.1 MM and Cross Validation
    2. 4.2 Double Kernel and Likelihood Cross Validation
    3. 4.3 Simulation Studies
  5. 5 Regularized and Monotone Quantile Regression
    1. 5.1 Adaptive Lasso for Heterogeneous Variable Selection
    2. 5.2 Monotone Spline for Quantile Monotonicity Preservation
  6. 6 Discussion
  7. A Additional Proofs
  8. B Additional Tables and Figures

##  1 Introduction

Predicting varied quantiles of the response variable is desirable in diverse
domains,111e.g., [3]; [39]; [47]; [26] as quantiles closer to the tails are of
more interest (compared to 0.5) in some applications. Various quantile
estimation approaches (e.g., [9]; [11]; [17]) have been proposed to model non-
constant correlation patterns between the dependent and independent variables
across different outcome quantiles, of which [28]’s method remains foremost.
Quantile regression introduced by [28] not only excels in efficiently
portraying a comprehensive picture of the regressor(s)-regressand’s dependency
structure, but also stays more robust against outliers than ordinary
regression of the mean.

Denote $n,p$ as the sample size, covariate dimension (including the intercept)
respectively and suppose we have observed a sample ($\bm{x_{i}},y_{i}),\text{
where }\bm{x_{i}}=(x_{i1},x_{i2},\dots,x_{ip})^{T}\text{ and }i=1,2,\dots,n$
(We usually let $x_{i1}s=1$ to denote intercept terms). For any given quantile
level $q\in(0,1)$ and any specified regression function
$f(\bm{X},\bm{\theta})$,222$Q_{q}(Y|\bm{X_{i}})=f(\bm{X_{i}},\bm{\theta}),i=1,2,\dots,n$,
where $\bm{\theta}=(\theta_{1},\theta_{2},\dots,\theta_{p})^{T}$ is the
parameter vector (corresponding to quantile level $q$) to be estimated. let
$f_{i}(\bm{\theta}):\mathbb{R}^{p}\rightarrow\mathbb{R}=f(\bm{x_{i}},\bm{\theta})\text{
and }r_{i}(\bm{\theta})=y_{i}-f_{i}(\bm{\theta})$333The $i^{\text{th}}$
residual., $i=1,2,\dots,n$. [28] defined the regression quantile at $q$ as the
$p$-dimensional vector $\bm{\hat{\theta}}$ that minimizes the asymmetric
$L_{1}$ loss function

$L(\bm{\theta})=\sum_{i=1}^{n}\rho_{q}(r_{i}(\bm{\theta})),$ (1.1)

where

$\rho_{q}(r_{i})=|r_{i}|[q\mathbbm{1}_{\\{r_{i}\geq
0\\}}+(1-q)\mathbbm{1}_{\\{r_{i}<0\\}}]=qr_{i}-r_{i}\mathbbm{1}_{\\{r_{i}<0\\}},\quad
i=1,2,\dots,n.$ (1.2)

Compared to traditional least-squares estimators, the minimizers of Equation
1.1 are much more robust to an array of non-normal error distributions while
performing similarly well for Gaussian linear models [28].

The prevailing most well-adopted quantile regression package in R is probably
still quantreg developed by [27], which estimates empirical regression
coefficients for specified quantile order(s) via a default simplex
algorithm[25] and an interior point method[36] triggering the Frisch-Newton
algorithm for huge data sets.444For data sets with at least several thousand
observations, specifying method=‘fn’ or method=‘pfn’ (Frisch-Newton Interior
Point Method with Preprocessing) may save computation time significantly. [35]
advanced in questing robust accelerated quantile regression methods [36] by
conducting simulation studies on four prominent approaches555The four
approaches are Alternating Direction Method of Multipliers (ADMM) algorithms,
Majorize-Minimize/Minorize-Maximize (MM) algorithms, Coordinate Descent (CD)
algorithms, and algorithms adopted by quantreg. including MM algorithms666As
proposed by [20] for non-regularized quantile regression and by [23] for
adaptive lasso-regularized quantile regression. and quantreg methods and
finding out that MM produced comparably accurate parameter estimates markedly
faster than the other three methods in non-regularized quantile and
composite777In composite quantile regression, we estimate coefficients for a
sequence of quantile orders simultaneously. quantile regression.

Indeed, MM offers much more than speeding up empirical quantile regression
coefficient estimation. As will be illustrated by this article, MM can be
conveniently extended to parametrically model multiple quantiles
simultaneously for nonparametric quantile regression problems with complicated
nonlinear correlation, which has not yet been discussed by current literature
to the best of our knowledge. After introducing MM algorithms and the
pioneering framework put forward by [20] applying MM algorithms to quantile
regression in Section 2, we start off with linear quantile regression in
Section 3, modeling the regression coefficients for different quantile orders
separately and simultaneously by various parametric bases transformations
including logistic and natural spline. Section 4 generalizes our MM algorithms
to “nonparametric quantile regression” and implements existing widely adopted
nonparametric counterpart methods to corroborate our MM algorithms’
effectiveness. Section 5 applies MM to quantile regression with special
requirements including variable selection and monotonicity preservation.
Finally, Section 6 summarizes our achievements and points toward potential
future work.

##  2 MM Algorithms

###  2.1 An Overview of the MM Algorithm

Although termed “algorithm”, MM, which stands for Majorize-Minimize for
minimization problems and Minorize-Maximize for maximization problems, is in
fact a strategy for creating an easy-to-minimize/maximize surrogate function
that majorizes/minorizes the objective function. Given an objective function
$L(\bm{\theta}):\mathbb{R}^{p}\rightarrow\mathbb{R}$ that is difficult to be
minimized directly and a current iteration’s parameter estimate
$\bm{\theta}^{(t)}$ ($t\in\mathbb{N}$), a surrogate function
$Q(\bm{\theta|\theta}^{(t)}):\mathbb{R}^{p}\rightarrow\mathbb{R}$ is called a
majorizer of $L(\bm{\theta})$ at $\bm{\theta}^{(t)}$ if it satisfies

$Q(\bm{\theta|\theta}^{(t)})\geq L(\bm{\theta})\text{ for all }\bm{\theta}$
(2.1)

and in particular

$Q(\bm{\theta}^{(t)}|\bm{\theta}^{(t)})=L(\bm{\theta}^{(t)}).$ (2.2)

If the value $\bm{\theta}^{(t+1)}$ at which $Q(\bm{\theta|\theta}^{(t)})$ is
minimized can be conveniently calculated, then $Q(\bm{\theta|\theta}^{(t)})$
yields a successful MM algorithm which drives the objective function
$L(\bm{\theta})$ further downhill each iteration thanks to the wonderful
descent property

$L(\bm{\theta}^{(t+1)})\leq Q(\bm{\theta}^{(t+1)}|\bm{\theta}^{(t)})\leq
Q(\bm{\theta}^{(t)}|\bm{\theta}^{(t)})=L(\bm{\theta}^{(t)}).$ (2.3)

To construct suitable majorizing/minorizing surrogate functions and thus
significantly simplify difficult optimization problems, MM typically exploits
convexity/concavity and miscellaneous inequalities [22] such as Jensen’s
Inequality (which brings about the well-known EM algorithm, a special case of
MM), Supporting Hyperplanes, Quadratic Upper Bound [2], the Arithmetic-
Geometric Mean Inequality, and the Cauchy-Schwartz Inequality. The MM
algorithm enjoys implementation simplicity and remarkable stability guaranteed
by Equation 2.3 while possessing complementary computation speed merits [22]
compared to the widely-adopted Newton-Raphson algorithm, and has hence gained
application attention in diversified fields.888e.g., [38]; [21]; [19] ;[23].

###  2.2 An MM Algorithm in Quantile Regression

For any given quantile level $q\in(0,1)$, the loss function $L(\bm{\theta})$
in Equation 1.1 is hard to directly minimize since it may accommodate multiple
minima and its components Equation 1.2 are non-differentiable when $r_{i}=0$
[20]. An MM algorithm may be constructed to bypass the non-differentiability
problem by setting ([20]; [22])

$Q(\bm{\theta}|\bm{\theta}^{(t)})=\sum_{i=1}^{n}\zeta_{q}(y_{i}-f_{i}(\bm{\theta})|y_{i}-f_{i}(\bm{\theta}^{(t)}))=\sum_{i=1}^{n}\zeta_{q}(r_{i}(\bm{\theta})|r_{i}(\bm{\theta}^{(t)})),$
(2.4)

where

$\zeta_{q}(r_{i}(\bm{\theta})|r_{i}(\bm{\theta}^{(t)}))=\frac{1}{4}\left\\{\frac{r_{i}^{2}(\bm{\theta})}{|r_{i}(\bm{\theta}^{(t)})|}+(4q-2)r_{i}(\bm{\theta})+|r_{i}(\bm{\theta}^{(t)})|\right\\}.$
(2.5)

It can be shown that the surrogate function in Equation 2.4 is a majorizer of
$L(\bm{\theta})$ at $\bm{\theta}^{(t)}$ (We leave the proof to Proposition
A.1. in Appendix A), yet this MM method suffers from the vexing drawback that
$Q(\bm{\theta}|\bm{\theta}^{(t)})$ is not defined when one
$r_{i}(\bm{\theta}^{(t)})=0$, a realistic concern that should be settled as
some $f_{i}(\bm{\theta}^{(t)})$, $i\in\\{1,2,\dots,n\\}\text{ may well
converge to }y_{i}\text{ as }t\rightarrow\infty.$999e.g., in the simplest case
when all $f_{i}(\bm{\theta}^{(t)})$, $i=1,2,...,n$ denote the sample mean
$\mu^{(t)}$, $\mu^{(t)}\text{ should }\rightarrow\text{ one }y_{i}\text{ as
}t\rightarrow\infty$ if $nq\not\in\mathbb{N}$ [20].

[20] addressed this issue by introducing a close-approximate function
$L_{\epsilon}(\bm{\theta})$101010We found out that a perturbation $\epsilon$
in the order of $10^{-4}$ or smaller usually generates good results.
($\epsilon>0$ denotes a perturbation) of the actual loss function Equation 1.1

$L_{\epsilon}(\bm{\theta})=\sum_{i=1}^{n}\rho_{q}^{\epsilon}(r_{i}(\bm{\theta}))=\sum_{i=1}^{n}\left\\{\rho_{q}(r_{i}(\bm{\theta}))-\frac{\epsilon}{2}\ln(\epsilon+|r_{i}(\bm{\theta})|)\right\\}$
(2.6)

and inventing a corresponding majorizer

$Q_{\epsilon}(\bm{\theta}|\bm{\theta}^{(t)})=\sum_{i=1}^{n}\zeta_{q}^{\epsilon}(r_{i}(\bm{\theta})|r_{i}(\bm{\theta}^{(t)}))=\frac{1}{4}\sum_{i=1}^{n}\left\\{\frac{r_{i}^{2}(\bm{\theta})}{\epsilon+|r_{i}(\bm{\theta}^{(t)})|}+(4q-2)r_{i}(\bm{\theta})+{\operatorname*{const}}_{i}\right\\},$
(2.7)

where ${\operatorname*{const}}_{i}$s are chosen to ensure that

$\zeta_{q}^{\epsilon}(r_{i}(\bm{\theta}^{(t)})|r_{i}(\bm{\theta}^{t}))=\frac{1}{4}\left\\{\frac{r_{i}^{2}(\bm{\theta}^{(t)})}{\epsilon+|r_{i}(\bm{\theta}^{(t)})|}+(4q-2)r_{i}(\bm{\theta}^{(t)})+{\operatorname*{const}}_{i}\right\\}=\rho_{q}^{\epsilon}(r_{i}(\bm{\theta}^{(t)})),\>i=1,2,\dots,n.$
(2.8)

By definition, $Q_{\epsilon}(\bm{\theta}|\bm{\theta}^{(t)})$ satisfies
Equation 2.2 as
$\zeta_{q}^{\epsilon}(r_{i}(\bm{\theta}^{(t)})|r_{i}(\bm{\theta}^{(t)}))=\rho_{q}^{\epsilon}(r_{i}(\bm{\theta}^{(t)}))$
for all $i\in\\{1,2,\dots,n\\}$. Given current residual value
$r^{(t)}=r(\bm{\theta}^{(t)})$ at iteration t, we let

$h(r)=\zeta_{q}^{\epsilon}(r|r^{(t)})-\rho_{q}^{\epsilon}(r)=\frac{1}{4}\left\\{\frac{r^{2}}{\epsilon+|r^{(t)}|}-2r+{\operatorname*{const}}+2\epsilon\cdot\ln(\epsilon+|r|)\right\\}+r\mathbbm{1}_{\\{r<0\\}}.$

Then

$\displaystyle\text{For }r\geq
0,\;h^{\prime}(r)=\frac{1}{2}\left(\frac{r}{\epsilon+|r^{(t)}|}-1+\frac{\epsilon}{\epsilon+r}\right)=\frac{1}{2}\cdot\frac{r(r-|r^{(t)}|)}{(\epsilon+|r^{(t)}|)(\epsilon+r)};$
(2.9) $\displaystyle\text{For
}r<0,\;h^{\prime}(r)=\frac{1}{2}\left(\frac{r}{\epsilon+|r^{(t)}|}+1-\frac{\epsilon}{\epsilon-r}\right)=-\frac{1}{2}\cdot\frac{r(r+|r^{(t)}|)}{(\epsilon+|r^{(t)}|)(\epsilon-r)}.$
(2.10)

Hence $h^{\prime}(r)=0\iff r=0\text{ or }r=\pm|r^{(t)}|\iff r=0\text{ or
}r=\pm r^{(t)}$ since the sets $\\{|r^{(t)}|,-|r^{(t)}|\\}$ and
$\\{r^{(t)},-r^{(t)}\\}$ are the same. We can see from Equations 2.9 and 2.10
that $h^{\prime}(r)\leq 0\text{ for
}r\in(-\infty,-|r^{(t)}|]\cup[0,|r^{(t)}|]\text{ and }\geq 0\text{ for
}r\in[-|r^{(t)}|,0]\cup[|r^{(t)}|,\infty),$ i.e., $h(r)$ is non-increasing for
$r\in(-\infty,-|r^{(t)}|]\cup[0,|r^{(t)}|]\text{ and non-decreasing for
}r\in[-|r^{(t)}|,0]\cup[|r^{(t)}|,\infty).$ Observe that $h(r)$ is an even
function and thus by Equation 2.8,

$h(-r^{(t)})=h(r^{(t)})=0\iff h(-|r^{(t)}|)=h(|r^{(t)}|)=0\Rightarrow h(r)\geq
0\text{ for all }r.\qquad$

Hence,

$\zeta_{q}^{\epsilon}(r|r^{(t)})\geq\rho_{q}^{\epsilon}(r)\text{ for all
}r\Rightarrow Q_{\epsilon}(\bm{\theta}|\bm{\theta}^{(t)})\geq
L_{\epsilon}(\bm{\theta})\text{ for all }\bm{\theta}.$

By the above reasoning, $Q_{\epsilon}(\bm{\theta}|\bm{\theta}^{(t)})$
satisfies Equation 2.1 as well and thus indeed majorizes
$L_{\epsilon}(\bm{\theta})$ at $\bm{\theta}^{(t)}$.

When $f_{i}(\bm{\theta})$s are linear,
$Q_{\epsilon}(\bm{\theta}|\bm{\theta}^{(t)})$ would be quadratic in
$\bm{\theta}$ and the MM algorithm can proceed in simplicity by explicitly
solving for $\bm{\theta}^{(t+1)}$ that minimizes
$Q_{\epsilon}(\bm{\theta}|\bm{\theta}^{(t)})$ and updating the current
parameter estimate $\bm{\theta}^{(t)}$ to $\bm{\theta}^{(t+1)}$. When
$f_{i}(\bm{\theta})$s are nonlinear, a Gauss-Newton approach may be adopted to
drive $Q_{\epsilon}(\bm{\theta}|\bm{\theta}^{(t)})$ downhill111111In which
case the descent property in Equation 2.3 still holds and guarantees the
algorithm’s stability. [20].

##  3 Linear Quantile Regression

###  3.1 Estimating Different Quantile Coefficients Separately

We first analyze the most extensively discussed and interpretation-wise
simplest linear case in quantile regression, a scenario under which the
$r_{i}(\bm{\theta})$s in Equation 2.7 can be straightforwardly written as
$\bm{r(\theta)=y-X\theta}$ ($\bm{y},\,\bm{X}$ contains observed sample values)
for any given quantile level $q\in(0,1)$, where

$\bm{r(\theta)}=\begin{pmatrix}r_{1}(\bm{\theta})\\\ r_{2}(\bm{\theta})\\\
\vdots\\\ r_{n}(\bm{\theta})\end{pmatrix},\;\bm{y}=\begin{pmatrix}y_{1}\\\
y_{2}\\\ \vdots\\\
y_{n}\end{pmatrix},\;\bm{X}=\begin{pmatrix}\bm{{x_{1}}}^{T}\\\
\bm{{x_{2}}}^{T}\\\ \vdots\\\
\bm{{x_{n}}}^{T}\end{pmatrix}=\begin{pmatrix}x_{11}&x_{12}&\dots&x_{1p}\\\
x_{21}&x_{22}&\dots&x_{2p}\\\ \vdots&\vdots&\ddots&\vdots\\\
x_{n1}&x_{n2}&\dots&x_{np}\end{pmatrix}$ (3.1)

since $f_{i}(\bm{\theta})=\bm{x_{i}}^{T}\bm{\theta}\text{ for }i=1,2,\dots,n.$
Typically, the first column of $\bm{X}$ is a column of $1s$ denoting the
intercept.

Given $q\in(0,1)$ and current parameter estimate $\bm{\theta}^{(t)}$, let

$\bm{W}^{(t)}=\begin{pmatrix}\frac{1}{\epsilon+|r_{1}(\bm{\theta}^{(t)})|}&&&\\\
&\frac{1}{\epsilon+|r_{2}(\bm{\theta}^{(t)})|}&&\\\ &&\ddots&\\\
&&&\frac{1}{\epsilon+|r_{n}(\bm{\theta}^{(t)})|}\end{pmatrix},\;\bm{c}_{n\times
1}=\begin{pmatrix}4q-2\\\ 4q-2\\\ \vdots\\\ 4q-2\end{pmatrix}.$ (3.2)

Then Equation 2.7 becomes

$\displaystyle Q_{\epsilon}(\bm{\theta}|\bm{\theta}^{(t)})=$
$\displaystyle\frac{1}{4}\left(\bm{r(\theta)}^{T}\bm{W}^{(t)}\bm{r(\theta)}+\bm{c}^{T}\bm{r(\theta)}\right)+\operatorname*{const}$
$\displaystyle=$
$\displaystyle\frac{1}{4}\left(\bm{(y-X\theta)}^{T}\bm{W}^{(t)}\bm{(y-X\theta)}+\bm{c}^{T}\bm{(y-X\theta)}\right)+\operatorname*{const}.$
(3.3)

Thus minimizing $Q_{\epsilon}(\bm{\theta}|\bm{\theta}^{(t)})$ is equivalent to
minimizing

$g(\bm{\theta}|\bm{\theta}^{(t)})=\bm{(y-X\theta)}^{T}\bm{W}^{(t)}\bm{(y-X\theta)}+\bm{c}^{T}\bm{(y-X\theta)},$
(3.4)

which can be done directly as $g(\bm{\theta}|\bm{\theta}^{(t)})$ is quadratic
in $\bm{\theta}$. Setting

$\frac{\partial}{\partial\bm{\theta}}g(\bm{\theta}|\bm{\theta}^{(t)})=2\bm{X}^{T}\bm{W}^{(t)}\bm{X\theta}-2\bm{X}^{T}\bm{W}^{(t)}\bm{y}-\bm{X}^{T}\bm{c}=\bm{0}\,\text{
yields}$
$\bm{\theta}^{(t+1)}=(\bm{X}^{T}\bm{W}^{(t)}\bm{X})^{-1}\left(\bm{X}^{T}\bm{W}^{(t)}\bm{y}+\frac{1}{2}\bm{X}^{T}\bm{c}\right).$
(3.5)

$(\bm{X}^{T}\bm{W}^{(t)}\bm{X})^{-1}$ exists for a full rank $\bm{X}_{n\times
p}$ ($n>p$), as

$\displaystyle\bm{W}^{(t)}\text{ is diagonal whose entries are all strictly
positive }$ $\displaystyle\Rightarrow$ $\displaystyle\text{ for all }p\times
1\text{ vector }\bm{v}\neq\bm{0}_{p\times 1},$
$\displaystyle\bm{v}^{T}\bm{X}^{T}\bm{W}^{(t)}\bm{Xv}=\bm{(Xv)}^{T}\bm{W}^{(t)}\bm{(Xv)}>0\text{
since }\bm{Xv}\neq\bm{0}_{n\times 1}\text{ for a full rank }\bm{X}$
$\displaystyle\Rightarrow$ $\displaystyle\bm{X}^{T}\bm{W}^{(t)}\bm{X}\text{ is
positive definite for a full rank }\bm{X}$ $\displaystyle\Rightarrow$
$\displaystyle\bm{X}^{T}\bm{W}^{(t)}\bm{X}\text{ is invertible for a full rank
}\bm{X}.$

Since for a full rank $\bm{X}$, the Hessian matrix

$\frac{\partial^{2}}{\partial\bm{\theta}\partial\bm{\theta}^{T}}g(\bm{\theta}|\bm{\theta}^{(t)})=2\bm{X}^{T}\bm{W}^{(t)}\bm{X}$
(3.6)

is positive definite for all $\bm{\theta}$ and in particular
$\bm{\theta}^{(t+1)}$, $\bm{\theta}^{(t+1)}$ given in Equation 3.5 indeed
minimizes Equations 3.4 and 3.1 and gives the next iteration’s parameter
estimate for the MM algorithm.

We implemented this MM algorithm on the “pollution” data set [31] in R’s
SMPracticals package and compared our results with the ones produced by the rq
function in the quantreg package [27]. The data set consists of 60
observations and conveys information about a mortality measure ($y$) with 15
explanatory variables, out of which we considered $x_{1}$ (precipitation, or
“prec”), $x_{9}$ (percentage of non-white, or “nonw”), $x_{10}$ (percentage of
white collar, or “wwdrk”), and $x_{14}$ ($SO_{2}$ potential, or “so”), as
these four variables have been demonstrated by [33] to possess either fairly
constant significant ($x_{1},x_{9}$) or manifest varying ($x_{10},x_{14}$)
effects across different quantile orders.

Quantile Level | $\bm{0.1}$ | $\bm{0.3}$ | $\bm{0.5}$ | $\bm{0.7}$ | $\bm{0.9}$
---|---|---|---|---|---
Intercept | MM | 6.93236001846033 | 6.81129448032327 | 6.82937502691576 | 6.78261292514527 | 6.75103326513239
rq | 6.93236000893714 | 6.81129447876662 | 6.82937503030063 | 6.78261292566128 | 6.75103326659995
Prec | MM | 0.00176228378387 | 0.00176159829763 | 0.00219734102618 | 0.00256007248951 | 0.00379190890437
rq | 0.00176228385391 | 0.00176159830142 | 0.00219734120251 | 0.00256007251365 | 0.00379190894470
Nonw | MM | 0.00307198793255 | 0.00378855193952 | 0.00338337914962 | 0.00317285210519 | 0.00310426173527
rq | 0.00307198786522 | 0.00378855195678 | 0.00338337911208 | 0.00317285194220 | 0.00310426167488
Wwdrk | MM | -0.00564678398232 | -0.00247626438379 | -0.00285944783005 | -0.00170900792443 | -0.00108160523893
rq | -0.00564678381819 | -0.00247626435265 | -0.00285944806563 | -0.00170900789901 | -0.00108160528760
So | MM | 0.00051426030202 | 0.00040409834617 | 0.00037884373159 | 0.00035640682915 | 0.00019472676333
rq | 0.00051426030253 | 0.00040409834526 | 0.00037884377530 | 0.00035640682973 | 0.00019472676437

Table 1: Comparison of selected quantile orders’ coefficient estimates from rq
and MM (with number of iterations $N=100000$ and perturbation
$\epsilon=10^{-10}$) for linear quantile regression on
$y,x_{1},x_{9},x_{10},x_{14}$ of the “pollution” data set.

Reassuringly, we arrived at extremely close regression coefficient estimates
for various quantile order choices. A comparison of the two methods’ estimated
coefficients for quantile orders $q=0.1,\,0.3,\,0.5,\,0.7,\,0.9\,$ is
presented in Table 1. This lends us considerable confidence in the
effectiveness of this aforementioned MM algorithm (which also turned out to be
fast), and we are equipped to extend our MM algorithm to model linear
regression coefficients for different quantile orders concurrently rather than
separately.

###  3.2 Estimating Different Quantile Coefficients Simultaneously

It is of tremendous interest as well as utmost importance to efficiently
discern informative signals from and discard unwanted noises in the empirical
sample quantile functions. To date, most relevant efforts either resort to
varied smoothing techniques on the empirical quantile functions (e.g.,[7];
[8]) or adjust nonparametric fitting techniques to quantile regression
(e.g.,[46]), which more or less suffer from interpretation simplicity and
prediction efficiency. [16] addressed these concerns by bringing forward a
parsimonious parametric approach that directly models the linear regression
coefficients as smooth functions of $q$,121212With the aid of the Newton-
Raphson algorithm to fit models. which succeeds in effectively pooling
information across quantile levels. The method was later built into R via the
qrcm package [15].

Rather than estimating the regression coefficients directly but only one $q$ a
time as in the previous subsection, [16] introduced a basis of
$h\in\mathbb{N}$131313$h$ is small and in most cases less than 10 known
functions of the quantile order –
$\bm{b}(q)=(b_{1}(q),b_{2}(q),\dots,b_{h}(q))^{T},$141414$\bm{b}(q)$ is
required to induce a monotonically non-decreasing $Q(q|\bm{x,\theta})$ for
some $\bm{\theta}$. wrote the regression coefficient as
$\bm{\beta}_{q}=\bm{A}_{p\times h}\bm{b}(q)$ for all $q$,151515This way,
$\bm{\beta}_{q}$’s are represented as linear combinations of basis functions
of $q$ and are no longer estimated separately. and estimated instead the
$p\times h$ parameter matrix

$\bm{A}=\begin{pmatrix}\theta_{11}&\theta_{12}&\dots&\theta_{1h}\\\
\theta_{21}&\theta_{22}&\dots&\theta_{2h}\\\ \vdots&\vdots&\ddots&\vdots\\\
\theta_{p1}&\theta_{p2}&\dots&\theta_{ph}\end{pmatrix}.$ (3.7)

Given $\bm{A},$ the residual vector can now be expressed as a function of $q$
as well:
$\bm{r}(q)=(r_{1}(q),r_{2}(q),\dots,r_{n}(q))^{T}=\bm{y}-\bm{XAb}(q)$ , where
$\bm{y},\bm{X}$ are as in Equation 3.1.

Our MM algorithm described in Section 3.1 can be conveniently generalized to
fit the model in this section by rewriting the parameter matrix
$\bm{A}_{p\times h}$ in Equation 3.7 as

$\bm{\theta}_{hp\times
1}=(\theta_{11},\theta_{21},\dots,\theta_{p1},\dots\dots,\theta_{1h},\theta_{2h},\dots,\theta_{ph})^{T}$
(3.8)

to form the parameter vector to be estimated and replace the $\bm{X}_{n\times
p}$ in Section 3.1 by new design matrices $\bm{D}(q)_{n\times hp}$ for any
$q\in(0,1)$:

$\displaystyle\bm{D}(q)_{n\times hp}$ $\displaystyle=\bm{b}(q)^{T}_{1\times
h}\otimes\bm{X}_{n\times p}$ (3.9)
$\displaystyle=\begin{pmatrix}x_{11}b_{1}(q)&x_{12}b_{1}(q)&\dots&x_{1p}b_{1}(q)&\dots\dots&x_{11}b_{h}(q)&x_{12}b_{h}(q)&\dots&x_{1p}b_{h}(q)\\\
x_{21}b_{1}(q)&x_{22}b_{1}(q)&\dots&x_{2p}b_{1}(q)&\dots\dots&x_{21}b_{h}(q)&x_{22}b_{h}(q)&\dots&x_{2p}b_{h}(q)\\\
\vdots&\vdots&\ddots&\vdots&\ddots&\vdots&\vdots&\ddots&\vdots\\\
\vdots&\vdots&\ddots&\vdots&\ddots&\vdots&\vdots&\ddots&\vdots\\\
x_{n1}b_{1}(q)&x_{n2}b_{1}(q)&\dots&x_{np}b_{1}(q)&\dots\dots&x_{n1}b_{h}(q)&x_{n2}b_{h}(q)&\dots&x_{np}b_{h}(q)\end{pmatrix}.$

Then
$\bm{r}(q,\bm{\theta})=(r_{1}(q,\bm{\theta}),r_{2}(q,\bm{\theta}),\dots,r_{n}(q,\bm{\theta}))^{T}=\bm{y}-\bm{XAb}(q)=\bm{y}-\bm{D}(q)\bm{\theta}$.
Given $k\in\mathbb{N}$ quantile orders $q_{1},q_{2},\dots,q_{k}$ that we wish
to simultaneous simultaneously model, where $k$ is typically set large (e.g.,
100, 1000) for better performances, the surrogate function with perturbation
becomes

$\displaystyle Q_{\epsilon}(\bm{\theta}|\bm{\theta}^{(t)})=$
$\displaystyle\sum_{a=1}^{k}\sum_{i=1}^{n}\zeta_{q_{a}}^{\epsilon}(r_{i}(q_{a},\bm{\theta})|r_{i}(q_{a},\bm{\theta}^{(t)}))$
$\displaystyle=$
$\displaystyle\frac{1}{4}\sum_{a=1}^{k}\sum_{i=1}^{n}\left\\{\frac{r_{i}^{2}(q_{a},\bm{\theta})}{\epsilon+|r_{i}(q_{a},\bm{\theta}^{(t)})|}+(4q_{a}-2)r_{i}(q_{a},\bm{\theta})+{\operatorname*{const}}_{ia}\right\\}.$
(3.10)

Given current parameter estimate $\bm{\theta}^{(t)},\text{ for all
}a\in\\{1,2,\dots,k\\}$, let

$(\bm{W}_{a})_{n\times
n}^{(t)}=\begin{pmatrix}\frac{1}{\epsilon+|r_{1}(q_{a},\bm{\theta}^{(t)})|}&&&\\\
&\frac{1}{\epsilon+|r_{2}(q_{a},\bm{\theta}^{(t)})|}&&\\\ &&\ddots&\\\
&&&\frac{1}{\epsilon+|r_{n}(q_{a},\bm{\theta}^{(t)})|}\end{pmatrix},\;(\bm{c}_{a})_{n\times
1}=\begin{pmatrix}4q_{a}-2\\\ 4q_{a}-2\\\ \vdots\\\ 4q_{a}-2\end{pmatrix}.$
(3.11)

Then Section 3.2 can be written as

$\displaystyle Q_{\epsilon}(\bm{\theta}|\bm{\theta}^{(t)})=$
$\displaystyle\frac{1}{4}\sum_{a=1}^{k}\left\\{\bm{r}(q_{a},\bm{\theta})^{T}\bm{W}_{a}^{(t)}\bm{r}(q_{a},\bm{\theta})+\bm{c}_{a}^{T}\bm{r}(q_{a},\bm{\theta})+{\operatorname*{const}}_{a}\right\\}$
$\displaystyle=$
$\displaystyle\frac{1}{4}\sum_{a=1}^{k}\left\\{\bm{(y-D}(q_{a})\bm{\theta)}^{T}\bm{W}_{a}^{(t)}\bm{(y-D}(q_{a})\bm{\theta)}+\bm{c}_{a}^{T}\bm{(y-D}(q_{a})\bm{\theta)}+{\operatorname*{const}}_{a}\right\\},$
(3.12)

which is again quadratic in $\bm{\theta}$ and thus can be solved directly for
$\bm{\theta}^{(t+1)}$. Let

$g(\bm{\theta}|\bm{\theta}^{(t)})=\sum_{a=1}^{k}\left\\{\bm{(y-D}(q_{a})\bm{\theta)}^{T}\bm{W}_{a}^{(t)}\bm{(y-D}(q_{a})\bm{\theta)}+\bm{c}_{a}^{T}\bm{(y-D}(q_{a})\bm{\theta)}\right\\}.$
(3.13)

Then minimizing $g(\bm{\theta}|\bm{\theta}^{(t)})$ is equivalent to minimizing
$Q_{\epsilon}(\bm{\theta}|\bm{\theta}^{(t)})$.

$\displaystyle\frac{\partial}{\partial\bm{\theta}}g(\bm{\theta}|\bm{\theta}^{(t)})\bigg{|}_{\bm{\theta}=\bm{\theta}^{(t+1)}}=\sum_{a=1}^{k}\left\\{2\bm{D}^{T}(q_{a})\bm{W}_{a}^{(t)}\bm{D}(q_{a})\bm{\theta}^{(t+1)}-2\bm{D}^{T}(q_{a})\bm{W}_{a}^{(t)}\bm{y}-\bm{D}^{T}(q_{a})\bm{c}_{a}\right\\}=\bm{0}$
$\displaystyle\qquad\Rightarrow\bm{\theta}^{(t+1)}=\left[\sum_{a=1}^{k}\bm{D}^{T}(q_{a})\bm{W}_{a}^{(t)}\bm{D}(q_{a})\right]^{-1}\left[\sum_{a=1}^{k}\left(\bm{D}^{T}(q_{a})\bm{W}_{a}^{(t)}\bm{y}+\frac{1}{2}\bm{D}^{T}(q_{a})\bm{c}_{a}\right)\right],$
(3.14)

which is well-defined for a full rank $\bm{X}_{n\times p}$ ($n>p$), a small
$h\in\mathbb{N}$, and a big $k\in\mathbb{N}$. To see this, we note that
$\sum_{a=1}^{k}\bm{D}^{T}(q_{a})_{hp\times n}(\bm{W}_{a})^{(t)}_{n\times
n}\bm{D}(q_{a})_{n\times hp}$ can be written as $\bm{D}^{T}_{hp\times
nk}\bm{W}_{nk\times nk}^{(t)}\bm{D}_{nk\times hp}$,

where $\displaystyle\bm{D}_{nk\times hp}=\begin{pmatrix}\bm{D}(q_{1})_{n\times
hp}\\\ \bm{D}(q_{2})_{n\times hp}\\\ \vdots\\\ \bm{D}(q_{k})_{n\times
hp}\end{pmatrix}=\begin{pmatrix}\bm{b}(q_{1})^{T}_{1\times
h}\otimes\bm{X}_{n\times p}\\\ \bm{b}(q_{2})^{T}_{1\times
h}\otimes\bm{X}_{n\times p}\\\ \vdots\\\ \bm{b}(q_{k})^{T}_{1\times
h}\otimes\bm{X}_{n\times p}\end{pmatrix}=\begin{pmatrix}\bm{b}(q_{1})^{T}\\\
\bm{b}(q_{2})^{T}\\\ \vdots\\\ \bm{b}(q_{k})^{T}\end{pmatrix}_{k\times
h}\otimes\bm{X}_{n\times p}$ with
$\displaystyle\operatorname*{rank}(\bm{D})=\operatorname*{rank}\left(\begin{pmatrix}\bm{b}(q_{1})^{T}\\\
\bm{b}(q_{2})^{T}\\\ \vdots\\\
\bm{b}(q_{k})^{T}\end{pmatrix}\right)\times\operatorname*{rank}(\bm{X})=h\times
p=hp\quad\Rightarrow\bm{D}_{nk\times hp}\text{ is full-rank},$ and
$\displaystyle\bm{W}_{nk\times
nk}^{(t)}=\begin{pmatrix}(\bm{W}_{1})^{(t)}_{n\times n}&&&\\\
&(\bm{W}_{2})^{(t)}_{n\times n}&&\\\ &&\ddots&\\\
&&&(\bm{W}_{k})^{(t)}_{n\times n}\end{pmatrix}\text{ is diagonal and positive
definite }.$ Hence,
$\displaystyle\sum_{a=1}^{k}\bm{D}^{T}(q_{a})\bm{W}_{a}^{(t)}\bm{D}(q_{a})=\bm{D}^{T}\bm{W}^{(t)}\bm{D}\text{
is positive definite}$ (3.15) $\displaystyle\Rightarrow\;\;$
$\displaystyle\sum_{a=1}^{k}\bm{D}^{T}(q_{a})\bm{W}_{a}^{(t)}\bm{D}(q_{a})=\bm{D}^{T}\bm{W}^{(t)}\bm{D}\text{
is invertible.}$

By Section 3.2, $\bm{\theta}^{(t+1)}$ indeed minimizes Equations 3.13 and 3.2
since

$\frac{\partial^{2}}{\partial\bm{\theta}\partial\bm{\theta}^{T}}g(\bm{\theta}|\bm{\theta}^{(t)})=2\sum_{a=1}^{k}\bm{D}^{T}(q_{a})\bm{W}_{a}^{(t)}\bm{D}(q_{a})\text{
is positive definite for all }\bm{\theta}\text{ and in particular
}\bm{\theta}^{(t+1)}.$

With no additional assumptions, the parametric logistic basis

$\bm{b}(q)=(1,\ln(q),\ln(1-q))^{T}$ (3.16)

has been suggested to perform well generally by [16] and was thus adopted by
us. Since there are situations, though not frequent, when linear combinations
of Equation 3.16 drastically deviate from the actual coefficient
functions,161616See Section 3.3 when the errors follow a shifted rescaled
$\operatorname*{Beta}(0.5,0.5)$ distribution for instance. we also resorted to
flexible nonparametric bases. Natural splines (cubic splines that restrict
both tails to be linear, abbreviated as “ns”), which encompass both
flexibility retained from ordinary cubic splines and stability achieved by
imposing the tail linearity restriction [34] and have thus been extensively
utilized in regression analysis, were considered. We followed the truncated
power basis form presented by [18] and defined

$\displaystyle\bm{b}(q)_{h\times
1}=(1,q,S_{1}(q),S_{2}(q),\dots,S_{h-2}(q))^{T},\text{ where for all
}l\in\\{1,2,\dots,h-2\\},$ $\displaystyle
S_{l}(q)=(q-q_{l})_{+}^{3}-\frac{(q-q_{h-1})_{+}^{3}(q_{h}-q_{l})}{q_{h}-q_{h-1}}+\frac{(q-q_{h})_{+}^{3}(q_{h-1}-q_{l})}{q_{h}-q_{h-1}}$
(3.17)

for ns bases with $h>2\>(h\in\mathbb{N})$ knots
$q_{1},q_{2},\dots,q_{h}\in(0,1)$. A variety of knot choices were considered,
some equally spaced and some with shorter tails.171717Sensible as the tails
are restricted to be linear.

We implemented this section’s MM algorithm on the “pollution” data set with
variables $y,x_{1},x_{9},x_{10},x_{14}$ again and carried out the
corresponding quantreg, qrcm [15] routines. The results validated our MM
algorithm’s appropriateness once more. Figure 1 and Figures B.1, B.2, B.3,
B.4, B.5, B.6 and B.7 depict the estimated linear regression coefficients
$\hat{\bm{\beta}}(q)_{p\times 1}$ as functions of the quantile order
$q\in(0,1)$, from which we can see that MM’s estimates capture efficiently rq
empirical estimates’ trends and might be slightly better than their qrcm
counterparts for all models’ five coefficients. Figure 2 presents estimated
response values $\hat{y}(q|\bm{x})=\hat{Q}_{q}(Y|\bm{x})$ at $1000$ quantile
levels $q\in(0,1)$, taking $\bm{x}=(x_{1},x_{9},x_{10},x_{14})$ as the
observed sample mean vector. For all models here, MM effectively smooths the
empirical functions produced by rq and agrees well to qrcm. We also made such
conditional $\hat{y}(q|\bm{x})$ quantile plots for each of the $n=60$ sample
observations and got similar corroborative results. Those plots generated
assuming the logistic basis $\bm{b}(q)=(1,\ln(q),\ln(1-q))^{T}$ and natural
spline basis with knots $0.2,0.4,0.6,0.8$ are displayed in Figures B.8 and
B.9.

Figure 1: Estimated regression coefficients $\hat{\bm{\beta}}(q)$ versus
quantile levels $q$ plot for our models with the logistic basis and 3 natural
spline bases. The orange lines are rq’s empirical regression coefficient
estimates at 1000 quantile levels. The blue and green lines denote estimated
$\hat{\bm{\beta}}(q)$ produced by qrcm and our MM algorithm described in this
section, respectively. Figure 2: Quantile plots
($\hat{Q}_{q}(Y|\bm{x})=\hat{y}(q|\bm{x})$ versus $q$, $q\in(0,1)$) for six
regression models’ predicted $\hat{y}$, given the independent variables
$\bm{x}_{p\times 1}=\sum_{i=1}^{n}\bm{x}_{i}/n\>(p=4+1=5,\>n=60)$. We have
adopted the logistic basis (A) and natural spline bases with knots 0.1, 0.3,
0.5, 0.7, 0.9 (B), 0.2, 0.4, 0.6, 0.8 (C), 0.1, 0.3, 0.7, 0.9 (D), 0.25, 0.5,
0.75 (E), and 0.1, 0.5, 0.9 (F).

###  3.3 Simulation Studies

We now carry out simulation studies to answer the following questions:

  1. 1.

Are our MM estimates close to those obtained from rq, qrcm?

  2. 2.

Does estimating quantile coefficients simultaneously (the idea in Section 3.2)
improve model fit (compared to modeling quantile coefficients separately as in
Section 3.1)?

  3. 3.

When modeling quantile coefficients as linear combinations of basis functions,
do parametric (the logistic basis – Equation 3.16) or nonparametric (ns bases
– Section 3.2) basis functions perform better?

We adopted the heterogeneous linear regression model

$Y_{i}=1+3X_{i}+(1+2X_{i})e_{i},\;i=1,2,\dots,n=\text{ sample size},$ (3.18)

where the regressors
$X_{i}\overset{\operatorname*{iid}}{\sim}\operatorname*{Unif}(0,1)$ are one-
dimensional uniform random variables independent of the $e_{i}$s, and the
errors $e_{i}$ are IID following an array of distributions:

  1. 1.

$e_{i}\sim N(0,1)$;

  2. 2.

$e_{i}\sim t(10)$ (heavy-tailed);

  3. 3.

$e_{i}=5(B_{i}-\frac{1}{2})$, where
$B_{i}\sim\operatorname*{Beta}(\alpha,\beta)$ with different representative
shape parameters $\alpha,\beta$ (short-tailed):

    1. (a)

$B_{i}\sim\operatorname*{Beta}(2,2)$;

    2. (b)

$B_{i}\sim\operatorname*{Beta}(\frac{1}{2},\frac{1}{2})$ (bathtub-shaped);

    3. (c)

$B_{i}\sim\operatorname*{Beta}(2,5)$ (asymmetric);

  4. 4.

$e_{i}=\ln(\epsilon_{i}),\text{ where }\epsilon_{i}\sim\exp(1)$ (asymmetric).

The above distributions were chosen to represent differing theoretical
quantile coefficient functions (Figure 3) and conditional quantile curves
(Figure 4), so as to answer our aforementioned three questions more
comprehensively.

Figure 3: Theoretical quantile coefficient functions for the six error
distributions – standard normal (A), $t(10)$ (B), standard log exponential
(C), and short-tailed beta distributions $\operatorname*{Beta}(2,2)$ (D),
$\operatorname*{Beta}(0.5,0.5)$ (E), $\operatorname*{Beta}(2,5)$ (F). Figure
4: Theoretical conditional quantile curves at
$q=0.05,0.1,0.25,0.5,0.75,0.9,0.95$ for the six error distributions – standard
normal (A), $t(10)$ (B), standard log exponential (C), and short-tailed beta
distributions $\operatorname*{Beta}(2,2)$ (D), $\operatorname*{Beta}(0.5,0.5)$
(E), $\operatorname*{Beta}(2,5)$ (F).

With the exact relationships between $Y_{i}$ and $X_{i}$, $e_{i}$ (Equation
3.18) and theoretical distributions of $X_{i}$ and $e_{i}$ for all
$i\in\\{1,\ldots,n\\}$, we can derive the theoretical conditional quantile
function

$\displaystyle Q_{q}(Y|x)=F_{y}^{-1}(q|x)$
$\displaystyle=1+3x+(1+2x)Q_{q}(e)\text{ since }Q_{q}(e|x)=Q_{q}(e)\text{ by
independence}$ $\displaystyle=[1+Q_{q}(e)]+[3+2Q_{q}(e)]x,\>q\in(0,1).$ (3.19)

In the special case where
$e_{i}\overset{\operatorname*{iid}}{\sim}\text{standard
logistic},Q_{q}(e)=\ln(q)-\ln(1-q)\text{ and hence
}Q_{q}(Y|x)=[1+\ln(q)-\ln(1-q)]+[3+2(\ln(q)-\ln(1-q))]x\Rightarrow$ the
theoretical coefficients are exact linear combinations of logistic basis
functions. Since the normal distribution, which is largely similar to t
distributions despite being less heavy-tailed, is very close to the logistic
distribution, we expect the logistic basis to work well for normal and $t(10)$
errors. This indeed turned out to be the case (Figure 5, Tables 2 and 3).

We calculated our deviation statistic – the Integrated Mean Squared Error
(IMSE)

$\operatorname*{IMSE}=\frac{1}{N}\sum_{l=1}^{N}\sum_{i=1}^{n}\left(Q_{q}(Y|x_{i})-\hat{Q}_{q}^{(l)}(Y|x_{i})\right)^{2}$
(3.20)

using $N$ random samples, each of size $n$, at quantile levels
$q=0.1,\>0.2,\>\dots,\>0.9$ for a variety of linear quantile regression
approaches. These methods include rq from quantreg, qMM (our MM algorithm
described in Section 3.1), iqr from qrcm, and qrcmMM (our MM algorithm
detailed in Section 3.2) with the logistic basis and natural spline bases.

$e\sim N(0,1)$
---
Quantile Level | $\bm{0.1}$ | $\bm{0.2}$ | $\bm{0.3}$ | $\bm{0.4}$ | $\bm{0.5}$ | $\bm{0.6}$ | $\bm{0.7}$ | $\bm{0.8}$ | $\bm{0.9}$
rq | 23.47573 | 15.02733 | 13.92952 | 12.05381 | 11.74653 | 13.63009 | 16.08107 | 18.41651 | 23.59899
qMM | 23.38479 | 15.02880 | 13.89654 | 12.03426 | 11.69038 | 13.63159 | 16.07505 | 18.41314 | 23.57543
qrcmlogistic | 16.63912 | 11.33287 | 10.32083 | 10.14422 | 10.27942 | 10.61187 | 11.20440 | 12.63035 | 19.28948
MMlogistic | 20.71841 | 11.74246 | 10.43766 | 10.42290 | 10.61639 | 10.80895 | 11.21838 | 13.10416 | 24.29688
qrcmspline 0.1_0.5_0.9 | 19.79880 | 14.36819 | 13.49741 | 11.70033 | 10.80515 | 11.58679 | 13.99664 | 16.10136 | 21.29121
MMspline 0.1_0.5_0.9 | 26.34934 | 12.91527 | 11.67682 | 11.05815 | 10.82807 | 11.14313 | 12.46402 | 14.82816 | 27.60851
qrcmspline 0.2_0.5_0.8 | 19.77772 | 14.40166 | 13.40390 | 11.76487 | 11.02821 | 11.62793 | 13.93497 | 16.19392 | 21.26600
MMspline 0.2_0.5_0.8 | 26.19061 | 12.96204 | 11.61984 | 11.14272 | 11.06685 | 11.20513 | 12.41433 | 14.86081 | 27.32316
qrcmspline 0.1_0.3_0.7_0.9 | 18.63382 | 12.20082 | 10.96115 | 11.17010 | 11.16893 | 12.71730 | 12.89961 | 14.15499 | 20.97778
MMspline 0.1_0.3_0.7_0.9 | 22.22681 | 12.80523 | 11.31631 | 11.13256 | 11.20563 | 12.46787 | 13.12265 | 14.55093 | 23.70929
qrcmspline 0.2_0.4_0.6_0.8 | 18.34301 | 13.39207 | 10.43786 | 11.85352 | 11.46597 | 13.74518 | 12.10325 | 15.33940 | 20.55781
MMspline 0.2_0.4_0.6_0.8 | 22.18026 | 13.18196 | 11.00228 | 11.60724 | 11.51331 | 13.18410 | 12.54544 | 14.89629 | 23.73108
qrcmspline 0.1_0.3_0.5_0.7_0.9 | 18.61639 | 12.79785 | 11.90895 | 11.01052 | 11.69196 | 12.58550 | 14.17558 | 15.16782 | 21.07862
MMspline 0.1_0.3_0.5_0.7_0.9 | 22.13956 | 13.35260 | 12.37101 | 10.87306 | 11.67903 | 12.38913 | 14.66483 | 15.39758 | 23.54844
$e\sim t(10)$
Quantile Level | $\bm{0.1}$ | $\bm{0.2}$ | $\bm{0.3}$ | $\bm{0.4}$ | $\bm{0.5}$ | $\bm{0.6}$ | $\bm{0.7}$ | $\bm{0.8}$ | $\bm{0.9}$
rq | 33.05197 | 20.54747 | 16.46722 | 14.99517 | 14.39148 | 13.68401 | 15.77772 | 18.63774 | 32.75290
qMM | 32.97132 | 20.45028 | 16.46175 | 15.00168 | 14.39526 | 13.70462 | 15.77588 | 18.60911 | 32.80304
qrcmlogistic | 26.71196 | 15.63156 | 12.77272 | 11.82905 | 11.49945 | 11.50780 | 12.10392 | 14.50295 | 24.53814
MMlogistic | 32.71580 | 16.66133 | 13.06830 | 12.15440 | 11.89980 | 11.86337 | 12.39448 | 15.27617 | 29.16942
qrcmspline 0.1_0.5_0.9 | 29.86629 | 21.29165 | 17.69877 | 13.45689 | 12.00307 | 13.91446 | 17.81342 | 20.16257 | 28.08354
MMspline 0.1_0.5_0.9 | 38.80241 | 18.57822 | 15.07723 | 12.81900 | 12.24541 | 13.16546 | 14.99494 | 17.18600 | 36.62121
qrcmspline 0.2_0.5_0.8 | 29.71517 | 21.44759 | 17.69246 | 13.49666 | 12.19712 | 13.96945 | 17.77214 | 20.26652 | 27.93254
MMspline 0.2_0.5_0.8 | 38.55690 | 18.68881 | 15.09070 | 12.89724 | 12.45744 | 13.21015 | 14.93143 | 17.24016 | 36.26492
qrcmspline 0.1_0.3_0.7_0.9 | 30.60296 | 18.54073 | 14.97582 | 14.32696 | 12.55400 | 13.48449 | 13.77033 | 17.43950 | 30.18305
MMspline 0.1_0.3_0.7_0.9 | 31.97637 | 18.19442 | 14.99380 | 13.72525 | 12.41510 | 13.09953 | 13.98424 | 16.65078 | 29.96764
qrcmspline 0.2_0.4_0.6_0.8 | 28.80200 | 20.55261 | 14.16197 | 15.60974 | 12.87046 | 14.55130 | 13.15664 | 19.72216 | 28.42172
MMspline 0.2_0.4_0.6_0.8 | 31.58774 | 18.61992 | 14.45499 | 14.57628 | 12.69555 | 13.78936 | 13.65927 | 17.45430 | 30.24770
qrcmspline 0.1_0.3_0.5_0.7_0.9 | 31.45977 | 18.64241 | 15.52735 | 14.40521 | 13.78792 | 13.57064 | 14.30894 | 17.47323 | 30.67790
MMspline 0.1_0.3_0.5_0.7_0.9 | 32.97841 | 18.19899 | 15.68273 | 13.93089 | 13.83919 | 13.14542 | 14.67902 | 17.06956 | 31.89505
$e=\ln(\epsilon),\textbf{ where }\epsilon\sim\exp(1)$
Quantile Level | $\bm{0.1}$ | $\bm{0.2}$ | $\bm{0.3}$ | $\bm{0.4}$ | $\bm{0.5}$ | $\bm{0.6}$ | $\bm{0.7}$ | $\bm{0.8}$ | $\bm{0.9}$
rq | 73.37499 | 38.26630 | 26.81722 | 17.52223 | 14.29766 | 12.30936 | 10.84083 | 10.59826 | 11.56988
qMM | 73.21025 | 38.25719 | 26.81489 | 17.51828 | 14.27827 | 12.32188 | 10.85397 | 10.59170 | 11.52108
qrcmlogistic | 25.43813 | 12.62038 | 8.62886 | 6.97765 | 6.14378 | 5.60514 | 5.15800 | 4.77733 | 5.14984
MMlogistic | 31.29570 | 13.86233 | 8.96239 | 7.21636 | 6.42851 | 5.85314 | 5.24740 | 4.83084 | 7.07683
qrcmspline 0.1_0.5_0.9 | 32.20710 | 23.84106 | 18.83181 | 9.73063 | 6.46465 | 7.09937 | 7.18537 | 5.49625 | 7.07096
MMspline 0.1_0.5_0.9 | 52.16195 | 16.56177 | 13.53939 | 8.72521 | 6.58719 | 6.31578 | 5.68999 | 4.91736 | 10.82732
qrcmspline 0.2_0.5_0.8 | 32.54457 | 24.79630 | 19.51990 | 9.63995 | 6.61376 | 7.20781 | 6.78265 | 5.23567 | 6.90289
MMspline 0.2_0.5_0.8 | 53.18192 | 16.89296 | 13.87613 | 8.67291 | 6.68429 | 6.36447 | 5.48888 | 4.89601 | 10.33078
qrcmspline 0.1_0.3_0.7_0.9 | 31.27132 | 20.49207 | 10.11653 | 9.64164 | 7.17527 | 6.34138 | 6.31523 | 4.38597 | 5.69953
MMspline 0.1_0.3_0.7_0.9 | 33.12593 | 16.50738 | 10.48543 | 8.62139 | 6.80724 | 6.07912 | 5.74202 | 4.65836 | 6.02686
qrcmspline 0.2_0.4_0.6_0.8 | 28.54116 | 25.45418 | 10.32736 | 10.74160 | 7.42589 | 7.17197 | 5.61831 | 4.47322 | 5.48634
MMspline 0.2_0.4_0.6_0.8 | 35.94247 | 17.94113 | 10.77236 | 9.08317 | 6.93095 | 6.51434 | 5.41351 | 4.64086 | 6.35627
qrcmspline 0.1_0.3_0.5_0.7_0.9 | 34.80144 | 18.07100 | 11.89301 | 9.57924 | 7.43537 | 6.68418 | 5.84666 | 5.03552 | 5.31413
MMspline 0.1_0.3_0.5_0.7_0.9 | 33.58544 | 16.07549 | 11.33470 | 8.63209 | 7.38014 | 6.14090 | 5.79588 | 4.90515 | 6.64794

Table 2: IMSE (as in Equation 3.20) calculated when the errors follow a
standard normal distribution (as in Item 1) with $N=300$, $n=100$, a heavy-
tailed $t$ distribution with df 10 (as in Item 2) and $N=650,n=100$, or an
asymmetric log exponential distribution (as in Item 4) with $N=300$, $n=200$.
The methods qrcmlogistic and MMlogistic denote respectively qrcm and qrcmMM
with the logistic basis Equation 3.16; the methods qrcmspline and MMspline
denote qrcm and qrcmMM with natural spline bases Section 3.2 containing the
specified knots. We used the test quantiles $q=0.1,0.2,\dots,0.9$ to fit all
MM models concerned.

$e=5(B-0.5),\textbf{ where }B\sim\operatorname*{Beta}(2,2)$
---
Quantile Level | $\bm{0.1}$ | $\bm{0.2}$ | $\bm{0.3}$ | $\bm{0.4}$ | $\bm{0.5}$ | $\bm{0.6}$ | $\bm{0.7}$ | $\bm{0.8}$ | $\bm{0.9}$
rq | 19.40513 | 21.24836 | 21.84782 | 20.25034 | 20.38719 | 19.86309 | 20.99954 | 18.94406 | 18.82771
qMM | 19.39963 | 21.20717 | 21.84134 | 20.24341 | 20.37637 | 19.85665 | 20.99925 | 18.90203 | 18.81788
qrcmlogistic | 17.28585 | 23.65278 | 24.25503 | 19.47401 | 16.08083 | 16.56463 | 19.47630 | 19.21028 | 15.50166
MMlogistic | 26.08059 | 18.42083 | 20.64423 | 18.55636 | 16.34528 | 16.03993 | 16.69002 | 15.17113 | 24.60922
qrcmspline 0.1_0.5_0.9 | 18.77248 | 15.04018 | 15.37462 | 15.70386 | 16.71546 | 17.66714 | 17.81084 | 15.82070 | 16.45175
MMspline 0.1_0.5_0.9 | 28.00586 | 15.48621 | 14.52098 | 15.69961 | 17.17998 | 17.24582 | 16.05263 | 15.16452 | 26.12272
qrcmspline 0.2_0.5_0.8 | 18.82177 | 15.01757 | 15.19645 | 15.83508 | 17.18561 | 17.83420 | 17.59456 | 15.72016 | 16.52196
MMspline 0.2_0.5_0.8 | 27.80989 | 15.53850 | 14.43131 | 15.86777 | 17.65450 | 17.37451 | 15.83232 | 15.12501 | 25.98825
qrcmspline 0.1_0.3_0.7_0.9 | 15.61443 | 14.26121 | 18.16157 | 18.79261 | 16.98778 | 17.67009 | 17.25622 | 14.15883 | 13.87907
MMspline 0.1_0.3_0.7_0.9 | 21.13871 | 15.45801 | 18.16402 | 18.52007 | 17.15530 | 17.70690 | 17.39621 | 14.83075 | 19.21747
qrcmspline 0.2_0.4_0.6_0.8 | 15.89827 | 14.07193 | 16.55244 | 20.68768 | 17.67756 | 19.07520 | 16.27214 | 14.45947 | 14.09082
MMspline 0.2_0.4_0.6_0.8 | 21.54050 | 14.99390 | 16.81811 | 19.85447 | 17.88920 | 18.86493 | 16.63567 | 14.71878 | 19.78943
qrcmspline 0.1_0.3_0.5_0.7_0.9 | 15.05536 | 15.73193 | 19.31529 | 19.01037 | 19.29097 | 17.84373 | 18.59047 | 15.47126 | 13.60727
MMspline 0.1_0.3_0.5_0.7_0.9 | 19.73116 | 16.46445 | 20.29608 | 18.87562 | 19.32885 | 17.96031 | 19.16456 | 16.01563 | 19.63744
$e=5(B-0.5),\textbf{ where }B\sim\operatorname*{Beta}(0.5,0.5)$
Quantile Level | $\bm{0.1}$ | $\bm{0.2}$ | $\bm{0.3}$ | $\bm{0.4}$ | $\bm{0.5}$ | $\bm{0.6}$ | $\bm{0.7}$ | $\bm{0.8}$ | $\bm{0.9}$
rq | 6.20090 | 30.16600 | 66.34303 | 107.05335 | 116.68954 | 99.82682 | 71.25583 | 32.91593 | 5.70905
qMM | 6.19859 | 30.14163 | 66.28704 | 107.03958 | 116.43953 | 99.84927 | 71.24012 | 32.99189 | 5.70775
qrcmlogistic | 62.42500 | 158.62156 | 252.99420 | 162.63501 | 98.99562 | 167.44816 | 258.69026 | 160.01989 | 66.08634
MMlogistic | 33.54415 | 191.80439 | 279.24975 | 171.89278 | 99.68156 | 173.93595 | 281.85207 | 192.79201 | 33.90016
qrcmspline 0.1_0.5_0.9 | 12.30912 | 56.46907 | 108.53049 | 110.44380 | 103.77000 | 112.89712 | 110.91832 | 56.91800 | 13.33990
MMspline 0.1_0.5_0.9 | 8.46361 | 50.35202 | 105.23387 | 110.84362 | 105.79604 | 113.44203 | 108.08538 | 51.48711 | 9.54160
qrcmspline 0.2_0.5_0.8 | 12.59758 | 55.20622 | 106.72386 | 112.13637 | 108.35890 | 114.49749 | 108.92259 | 55.58130 | 13.54125
MMspline 0.2_0.5_0.8 | 8.59840 | 48.84185 | 103.30063 | 112.71954 | 110.70629 | 115.08204 | 105.77929 | 49.80030 | 9.73125
qrcmspline 0.1_0.3_0.7_0.9 | 8.69627 | 32.37055 | 61.21867 | 75.60995 | 76.51542 | 75.89446 | 61.90748 | 32.90727 | 8.56311
MMspline 0.1_0.3_0.7_0.9 | 6.92579 | 31.66823 | 65.79468 | 82.51446 | 83.84504 | 81.18353 | 64.89839 | 31.98491 | 7.32607
qrcmspline 0.2_0.4_0.6_0.8 | 8.83155 | 35.91559 | 62.05128 | 94.73341 | 85.69125 | 92.82651 | 61.72181 | 36.25315 | 8.81373
MMspline 0.2_0.4_0.6_0.8 | 6.45397 | 32.33890 | 64.27563 | 95.64082 | 91.11935 | 92.61981 | 62.93832 | 32.44364 | 6.97216
qrcmspline 0.1_0.3_0.5_0.7_0.9 | 8.82865 | 31.98735 | 58.10374 | 85.39182 | 106.84981 | 84.26495 | 58.40157 | 32.94730 | 8.59773
MMspline 0.1_0.3_0.5_0.7_0.9 | 6.10987 | 32.21058 | 63.22233 | 88.90208 | 109.75912 | 87.29065 | 63.20782 | 32.78600 | 5.76422
$e=5(B-0.5),\textbf{ where }B\sim\operatorname*{Beta}(2,5)$
Quantile Level | $\bm{0.1}$ | $\bm{0.2}$ | $\bm{0.3}$ | $\bm{0.4}$ | $\bm{0.5}$ | $\bm{0.6}$ | $\bm{0.7}$ | $\bm{0.8}$ | $\bm{0.9}$
rq | 5.34402 | 6.30304 | 6.79065 | 8.39454 | 9.94945 | 11.28485 | 13.64479 | 14.73701 | 22.63101
qMM | 5.34064 | 6.30523 | 6.78039 | 8.40128 | 9.94809 | 11.28231 | 13.63609 | 14.64570 | 22.74989
qrcmlogistic | 4.12343 | 6.45969 | 7.36264 | 7.41939 | 7.73562 | 8.70354 | 10.12316 | 11.46949 | 17.48520
MMlogistic | 5.88103 | 5.41413 | 6.95313 | 7.67399 | 8.14397 | 8.74280 | 9.58557 | 11.31847 | 24.60911
qrcmspline 0.1_0.5_0.9 | 5.06382 | 4.72333 | 6.87298 | 8.07839 | 8.13547 | 9.57448 | 14.54964 | 17.87481 | 18.52272
MMspline 0.1_0.5_0.9 | 7.69693 | 4.71670 | 6.05908 | 7.65647 | 8.47003 | 9.49322 | 11.91206 | 13.53661 | 27.80076
qrcmspline 0.2_0.5_0.8 | 4.97562 | 4.47069 | 6.52237 | 8.22355 | 8.41503 | 9.55688 | 14.94913 | 18.52643 | 18.61279
MMspline 0.2_0.5_0.8 | 7.37637 | 4.67351 | 5.88621 | 7.75560 | 8.67979 | 9.48549 | 12.06263 | 13.69319 | 28.13956
qrcmspline 0.1_0.3_0.7_0.9 | 3.81142 | 4.31881 | 6.96217 | 7.80176 | 8.86790 | 10.62212 | 10.44315 | 14.75794 | 20.45766
MMspline 0.1_0.3_0.7_0.9 | 5.57674 | 4.73463 | 6.71543 | 7.83245 | 8.94939 | 10.60128 | 11.18215 | 13.21737 | 20.94692
qrcmspline 0.2_0.4_0.6_0.8 | 3.91619 | 4.00382 | 6.42124 | 8.47827 | 9.33310 | 11.70857 | 10.24959 | 17.66310 | 18.99403
MMspline 0.2_0.4_0.6_0.8 | 5.62832 | 4.55882 | 6.29610 | 8.24921 | 9.26629 | 11.21357 | 10.87005 | 13.73473 | 20.90892
qrcmspline 0.1_0.3_0.5_0.7_0.9 | 3.69995 | 4.60775 | 6.44292 | 8.09281 | 9.65834 | 10.67241 | 11.67749 | 13.58901 | 21.88336
MMspline 0.1_0.3_0.5_0.7_0.9 | 5.72442 | 5.02931 | 6.70284 | 7.85597 | 9.45304 | 10.43555 | 12.21406 | 13.20401 | 22.04415

Table 3: IMSE (as in Equation 3.20) calculated when the errors are distributed
as in Item 3 (rescaled shifted Beta distributions with different
representative shape parameters), with number of samples $N=300$ and sample
sizes $n=200$ for all three cases. The methods qrcmlogistic and MMlogistic
denote respectively qrcm and qrcmMM with the logistic basis Equation 3.16; the
methods qrcmspline and MMspline denote qrcm and qrcmMM with natural spline
bases Section 3.2 containing the specified knots. We used the test quantiles
$q=0.1,0.2,\dots,0.9$ to fit all MM models concerned. Figure 5: Mean estimated
quantile coefficient functions (dashed lines) using the logistic basis for the
six error distributions – standard normal (A), $t(10)$ (B), standard log
exponential (C), and short-tailed beta distributions
$\operatorname*{Beta}(2,2)$ (D), $\operatorname*{Beta}(0.5,0.5)$ (E),
$\operatorname*{Beta}(2,5)$ (F). The solid lines are the theoretical quantile
coefficient functions.

As shown in Tables 2 and 3, qMM IMSEs correspond extremely closely to their rq
counterparts at all tested quantile levels for all six aforementioned $e_{i}$
distributions. qrcmMM and qrcm are also comparable in general, especially at
non-extreme quantiles.181818Some bigger deviations at the tails may be due to
different model fitting quantile levels. We only used the 9 test quantiles
$q=0.1,0.2,\dots,0.9$ to fit all our MM models, whereas [16] typically used
250 quantile levels. Hence, we now focus on analyzing MM’s results. As
expected, qrcm with the logistic basis significantly outperforms qMM as well
as qrcm with nonparametric bases for normal, $t(10)$, and log exponential
errors. For the short-tailed $\operatorname*{Beta}(0.5,0.5)$ error
distribution, however, qrcm with the logistic basis loses out. As Figure 5E
depicts, qrcm with the logistic basis is indeed inadequate for short-tailed
$\operatorname*{Beta}(0.5,0.5)$ errors.191919Probably because linear
combinations of functions in the logistic basis concave up whereas the
theoretical quantile coefficient functions for short-tailed
$\operatorname*{Beta}(0.5,0.5)$ errors concave down. Under such scenarios when
the parametric logistic basis severely misspecifies the actual correlation
patterns, flexible nonparametric ns bases with more than three knots closely
estimate theoretical quantile coefficient functions (Figure 6) and do clearly
improve upon separate quantile estimates (Table 3), gaining us assurance in
the effectiveness of modeling quantile coefficients as linear combinations of
basis functions of $q$.

Figure 6: Mean estimated quantile coefficient functions (dashed lines) for the
short-tailed $\operatorname*{Beta}(0.5,0.5)$ error distribution. A: using ns
basis with knots $0.1,0.3,0.7,0.9$; B: using ns basis with knots
$0.1,0.3,0.5,0.7,0.9$. The solid lines are the theoretical quantile
coefficient functions.

##  4 Nonparametric Quantile Regression

Despite its popularity and general competence, the typical linear quantile
regression is inappropriate for data sets exhibiting manifest complex
nonlinear patterns. Quantile regression for universal nonlinear data, however,
has rarely been discussed. The handful of existing approaches in current
literature (e.g., [46]) rely on broader nonparametric methods like Kernel
Estimation adapted to quantile regression instead of parametric methods
convenient for interpretation and prediction. This section attempts to bridge
this gap by extending our MM algorithm constructed in Section 3.2 to
parametrically model this type of quantile regression problem, which quantreg
and qrcm algorithms are unable to handle. We start with two real bivariate
data sets displaying apparent complicated nonlinear correlations (Figure 7).

Figure 7: Scatter plots for two bivariate data sets. The “triceps” data set
(A) documents age and triceps skinfold thickness of 892 females under the age
of 50 in West Africa’s three villages. We used $\ln(\text{triceps})$ as the
response variable and took exponential of the estimated values to get back
triceps thickness’s original scale. The “WTD” data set (B) records 307 areas’
Water Table Depth (WTD) and rainfall flux at one site in South Africa’s Cape
Town. We transformed the variables as described by [32], i.e., take the
absolute value of WTD and let Rescaled Flux = $\frac{1}{7.3}(\ln(\text{Flux
}+6.9)+2.4)$.

###  4.1 MM and Cross Validation

It is in fact pretty straightforward to achieve the desired generalization via
MM. The first step is fairly customary in ordinary regression settings –
replace the original sample observation matrix $\bm{X}_{n\times p}$ ($p=1+1=2$
here for the two bivariate data sets) as in Equation 3.1 by some
transformations of the original variables, i.e.,

$\bm{X}=(\bm{1}_{n},\bm{x})=\begin{pmatrix}1&x_{1}\\\ 1&x_{2}\\\
\vdots&\vdots\\\
1&x_{n}\end{pmatrix}\overset{f_{1},f_{2},\dots,f_{p^{\prime}-1}}{\longrightarrow}\bm{X}_{new}=\begin{pmatrix}1&f_{1}(x_{1})&f_{2}(x_{1})&\dots&f_{p^{\prime}-1}(x_{1})\\\
1&f_{1}(x_{2})&f_{2}(x_{2})&\dots&f_{p^{\prime}-1}(x_{2})\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\
1&f_{1}(x_{n})&f_{2}(x_{n})&\dots&f_{p^{\prime}-1}(x_{n})\end{pmatrix}_{n\times
p^{\prime}},$ (4.1)

where $f_{1},f_{2},\dots,f_{p^{\prime}-1}$ are a set of chosen functions. With
this new design matrix $\bm{X}_{new}$, we can proceed to calculate the
corresponding $\bm{D}_{new}(q)_{n\times hp^{\prime}}$s in Equation 3.9 and
carry out our MM algorithm as delineated in Section 3.2.

As the correlation types are not specified and the two scatter plots displayed
in Figure 7 do not suggest any particular simple nonlinear patterns, we would
need highly flexible transformations applicable to universal nonlinear data
sets. We hence considered the flexible, stable, and extensively utilized
natural spline transformations [34] again, which were applied to covariates
instead of quantile orders this time. Since we first flexibly transform our
covariate vector $\bm{x}$ via the nonparametric natural spline transformation,
we name the subject for this section “nonparametric quantile regression”.

We conducted our MM algorithms on the two data sets with diverse combinations
of explanatory variable knots and quantile order bases. MM did appear to
smooth out rq empirical estimates reasonably well in all quantile plots of
$\hat{y}(q|x)=\hat{Q}_{q}(Y|x)$ versus $q$ that take $x$ as the sample mean
(Figures 9 and 11). Since the two data sets are bivariate, we can also
efficiently visualize MM’s prediction accuracy and smoothing effects via
estimated conditional quantile curves. As shown in Figures B.10 and 10 and
Figures B.11, B.12, 8, B.13, B.14, B.15, B.16, B.17, B.18 and B.19 plotting
$\hat{y}(x|q)$ against the concentrated $x$ values at selected representative
quantile levels $q=0.05,0.1,0.25,0.5,0.75,0.9,0.95$, our MM algorithm did a
good job in efficiently capturing the scatters’ trends and smoothing the
empirical estimates produced by rq for all combination choices.

“triceps” | “WTD”
---|---
$x$ $q$ | log | ns seq4 | ns seq7 |  $x$ $q$ | log | ns seq4 | ns asym4
seq3 | 760.265 | 761.903 | 761.582 | seq3 | 9313.770 | 9354.802 | 9350.683
seq4 | 744.874 | 746.292 | 746.893 | seq4 | 9624.973 | 9587.168 | 9588.477
seq5 | 747.905 | 748.785 | 748.992 | seq5 | 9431.753 | 9488.172 | 9556.216
seq6 | 747.635 | 748.080 | 748.352 | seq6 | 10446.993 | 10533.425 | 10522.805
seq7 | 741.567 | 742.729 | 743.431 | seq7 | 13292.171 | 13123.963 | 13037.342
seq8 | 739.630 | 740.801 | 741.479 | asym7 | 9354.379 | 9349.327 | 9351.190

Table 4: CV loss (Equation 4.2) of the “triceps” and “WTD” data sets for
selected MM $x,q$ settings. For natural spline transformations on $x$,
“seq$h$” $(h\in\mathbb{N})$ denotes putting $h$ equally spaced knots between
the minimum and maximum values of observed $x$’s and “asym7” for “WTD”
represents the knots $(0.2,0.4,0.5,0.55,0.6,0.7,0.9)$. For quantile order $q$
bases, “log”, “ns seq4”, “ns seq7”, “ns asym4” respectively denote the
logistic basis and natural spline bases with 4/7 equally spaced knots and the
knots $(0.3,0.4,0.6,0.8)$. We used 99 quantile levels $0.01,0.02,\dots,0.99$
and $N=100$ iterations to fit each qrcmMMxspline model in this chart.

To gain a better picture of the relevant out-of-sample prediction powers for
different covariate knots and quantile bases choices, we performed 10-fold
cross-validation (CV) by randomly splitting the data sets into 10
approximately equal-sized subgroups
$G_{1},G_{2},\dots,G_{10}\subset\\{1,2,\dots,n\\}$, picking $J=9$ validation
quantiles $q_{1}=0.1,q_{2}=0.2,\ldots,q_{J}=0.9$, and calculating the loss

$l_{CV}=\sum_{g=1}^{10}\sum_{i\in
G_{g}}\sum_{j=1}^{J}\rho_{q_{j}}\left(y_{i}-\bm{x}_{i}^{T}\hat{\bm{\beta}}^{(-G_{g})}(q_{j})\right),$
(4.2)

where each $\bm{x}_{i}$ denotes the $i^{\text{th}}$ row of $\bm{X}_{new}$ in
Equation 4.1 and each $\hat{\bm{\beta}}^{(-G_{g})}(q_{j})$ represents the
estimated coefficient vector based on the other 9 subgroups for spline-
transformed covariates at quantile level $q_{j}$. The smaller $l_{CV}$, the
better the model. As Table 4 suggests, basis choice for $q$ does not seem to
matter as much as ns knots choice to transform $x$ does.

Figure 8: Fitted quantile curves for the “triceps” data set - natural spline
transformation on $x$ with 6 equally spaced knots. We have adopted the
logistic basis (A) and natural spline bases with 3 (B), 4 (C), 5 (D), and 6
(E) equally spaced knots for quantile levels. Figure 9: Quantile plots for the
“triceps” data set taking $x$ as the sample mean. We natural-spline-
transformed $x$ using 3 (A), 4 (B), 5 (C), 6 (D), 7 (E), and 8 (F) equally
spaced knots between 0 and 1. Corresponding to each covariate transformation,
we adopted the logistic basis and natural spline bases with 3, 4, 5, and 6
equally spaced knots for quantiles on top of rq. Figure 10: Fitted quantile
curves for the “WTD” data set - natural spline transformation on $x$ with 7
knots $0.2,0.4,0.5,0.55,0.6,0.7,0.9$. We have adopted the logistic basis (A)
and natural spline bases with 3 (B), 4 (C), 5 (D), and 6 (E) equally spaced
knots for quantile levels. Figure 11: Quantile plots for the “WTD” data set
taking $x$ as the sample mean. We natural-spline-transformed $x$ using 3 (A),
4 (B), 5 (C), 6 (D), 8 (F) equally spaced knots between 0 and 1 and the
asymmetric knots $(0.2,0.4,0.5,0.55,0.6,0.7,0.9)$ (E). Corresponding to each
covariate transformation, we adopted the logistic basis and natural spline
bases with 3, 4, 5, and 6 equally spaced knots for quantiles on top of rq.

###  4.2 Double Kernel and Likelihood Cross Validation

To see how well our MM algorithms perform compared to existing nonparametric
methods, we adopted widely-used Kernel Estimation techniques on the two
bivariate real data sets to predict $\hat{y}(x|q)$ at
$q=0.05,0.1,0.25,0.5,0.75,0.9,0.95.$ Given a finite sample with observed
points $(x_{i},y_{i})$, $i=1,2,\dots,n$ and a new point with known $x$, the
estimated conditional CDF of this point’s $y$ coordinate is of the form

$\displaystyle\hat{F}(y|x)$ $\displaystyle=\sum_{i=1}^{n}w_{i}(x,h_{1})\cdot
I(y_{i}\leq y),\text{ where}$ (4.3) $\displaystyle w_{i}(x,h_{1})$
$\displaystyle=\frac{K\left(\frac{x_{i}-x}{h_{1}}\right)}{\sum_{j=1}^{n}K\left(\frac{x_{j}-x}{h_{1}}\right)}.$

To smooth our estimation in the $y$ direction, we replaced the deterministic
form $I(y_{i}\leq y)$ in Equation 4.3 by a standard normal CDF
$N(y_{i},h_{2})=\Phi\left(\frac{y-y_{i}}{h_{2}}\right)$ and formed a double
kernel conditional CDF estimate

$\hat{\hat{F}}(y|x)=\sum_{i=1}^{n}w_{i}(x,h_{1})\cdot\Phi\left(\frac{y-y_{i}}{h_{2}}\right)=\sum_{i=1}^{n}\left[\frac{K\left(\frac{x_{i}-x}{h_{1}}\right)}{\sum_{j=1}^{n}K\left(\frac{x_{j}-x}{h_{1}}\right)}\cdot\Phi\left(\frac{y-y_{i}}{h_{2}}\right)\right].$
(4.4)

The estimated conditional quantile function is thus

$\hat{\hat{Q}}_{q}(y|x)=\hat{\hat{F}}^{-1}(q|x),\>q\in(0,1).$ (4.5)

We adopted the normal kernel $\phi$ as the kernel function $K$ in Equation 4.4
in light of the normal kernel’s various favorable mathematical properties and
proceeded to produce quantile curves of $\hat{y}(x|q)$ against $x$ at each
$q\in\\{0.05,0.1,0.25,0.5,0.75,0.9,0.95\\}$ after specifying a set of sensible
choices for the bandwidths $h_{1},h_{2}$. We discovered and confirmed by
Likelihood cross-validation (LCV) that different $h_{2}$ values generated
estimates that are extremely alike, if not identical. As such, we present
quantile curves with a set of $h_{1}$ choices and a fixed adequate
$h_{2}=10^{-4}$ in Figures 12 and 13. The figures suggest that $h_{1}=10,12.5$
for “triceps” and $h_{1}=0.15,0.2$ for “WTD” most effectively and smoothly
capture the scatters’ trends.

To decide optimal bandwidth choices for Double Kernel Estimation, we carried
out leave-one-out likelihood cross-validation to quantitatively judge
different bandwidths’ out-of-sample prediction capabilities. As the estimated
conditional PDF is given by

$\hat{\hat{f}}(y|x)=\sum_{i=1}^{n}w_{i}(x,h_{1})\cdot\phi\left(\frac{y-y_{i}}{h_{2}}\right),$
(4.6)

a leave-one-out LCV likelihood can be expressed as

$\mathcal{L}_{LCV}=\prod_{j=1}^{n}\hat{\hat{f}}^{(-j)}(y_{j}|x_{j}),$ (4.7)

where $\hat{\hat{f}}^{(-j)}(y_{j}|x_{j})$ denotes the estimated PDF for a new
observation with $x=x_{j}$ derived using the entire sample except the point
$(x_{j},y_{j})$. The larger $\mathcal{L}_{LCV}$, the better the method. As
Table 5 suggests, $h_{1}\text{ approximately }\in(9,13)$ for “triceps”, whose
sample $x\in(0,52)$, and $h_{1}\text{ approximately }\in(0.12,0.22)$ for
“WTD”, whose sample $x\in(0,1)$, would be effective $x$ kernel bandwidth
choices, corresponding well to what the quantile curves in Figures 12 and 13
indicate. The values of $h_{2}$, on the other hand, do not seem to affect
Double Kernel Estimation’s prediction power much as long as $h_{2}$ is
appropriately small.

“triceps” | “WTD”
---|---
$h_{1}$ $h_{2}$ | $10^{-4}$ | $10^{-5}$ | $10^{-6}$ |  $h_{1}$ $h_{2}$ | $10^{-3}$ | $10^{-4}$ | $10^{-5}$
0.5 | 6.04978E-277 | 6.04978E-277 | 6.04978E-277 | 0.01 | 6.41081E-223 | 6.41081E-223 | 6.41081E-223
2.5 | 3.44698E+53 | 3.44698E+53 | 3.44698E+53 | 0.05 | 2.19999E-102 | 2.19999E-102 | 2.19999E-102
5 | 1.48980E+111 | 1.48980E+111 | 1.48980E+111 | 0.1 | 3.51314E-75 | 3.51314E-75 | 3.51314E-75
7.5 | 1.67443E+130 | 1.67443E+130 | 1.67443E+130 | 0.15 | 9.16157E-73 | 9.16157E-73 | 9.16157E-73
10 | 6.95286E+136 | 6.95286E+136 | 6.95286E+136 | 0.2 | 1.14691E-74 | 1.14691E-74 | 1.14691E-74
12.5 | 8.20678E+135 | 8.20678E+135 | 8.20678E+135 | 0.25 | 1.20404E-76 | 1.20404E-76 | 1.20404E-76
15 | 3.26345E+130 | 3.26345E+130 | 3.26345E+130 | 0.3 | 3.63631E-78 | 3.63631E-78 | 3.63631E-78

Table 5: Scaled LCV likelihood (Equation 4.7) calculated using some suitable
bandwidth choices for the “triceps” and “WTD” data sets. We have scaled the
likelihood product in Equation 4.7 to make sure that all resulting values can
be expressed in R. Hence, it is the relative magnitude rather than the
absolute magnitude that counts here. Figure 12: Quantile curves for the
“triceps” data set estimated by Double Kernel with $h_{1}=2.5$ (A), $5$ (B),
$7.5$ (C), $10$ (D), $12.5$ (E), $15$ (F) and $h_{2}=10^{-4}$. Figure 13:
Quantile curves for the “WTD” data set estimated by Double Kernel with
$h_{1}=0.05$ (A), $0.1$ (B), $0.15$ (C), $0.2$ (D), $0.25$ (E), $0.3$ (F) and
$h_{2}=10^{-4}.$

###  4.3 Simulation Studies

The quantile curves produced by Double Kernel seem to appear “rougher” than
those obtained from MM. To more convincingly and objectively compare the two
approaches, we resort to simulation studies again. We considered the three-
parameter generalized gamma (GG) distribution, as it incorporates many common
fewer-parameter distributions (e.g., Gamma, Weibull, Exponential
distributions).

A generalized gamma random variable Y with parameters $\theta,\beta,k$, or
$\mu,\sigma,k$, where

$\begin{cases}\mu=\ln(\theta)+\frac{\ln(k)}{\beta},\qquad\sigma=\frac{1}{\beta\sqrt{k}}\\\
\theta=\frac{e^{\mu}}{k^{\sigma\sqrt{k}}},\qquad\qquad\beta=\frac{1}{\sigma\sqrt{k}}\end{cases},$
(4.8)

has probability density function

$f(y;\theta,\beta,k)=\frac{\beta}{\theta^{k\beta}\Gamma(k)}y^{k\beta-1}\exp\left\\{-\left(\frac{y}{\theta}\right)^{\beta}\right\\},\qquad
y>0$ (4.9)

and quantile function

$Q(q;\mu,\sigma,k)=e^{\mu}\left\\{\frac{r(q;k)}{k}\right\\}^{\sigma\sqrt{k}},\qquad
q\in(0,1),$ (4.10)

where $r(q;k)$ is the quantile function of the
$\operatorname*{Gamma}(\theta=1,k)$, or
$\operatorname*{GG}(\theta=1,\beta=1,k)$ distribution. It is virtually the
$\frac{1}{\beta}^{\text{th}}$ power of a
$\operatorname*{Gamma}\left(\theta^{\beta},k\right)$ random variable and can
be generated either by substituting a $\operatorname*{Unif}(0,1)$ random
variable as $q$ into Equation 4.10 or by generating a
$\operatorname*{Gamma}\left(\theta^{\beta},k\right)$ random variable and
taking the $\frac{1}{\beta}^{\text{th}}$ power of it. For implementation in
software like R, the first method is preferred, as the second may fail for
some parameter choices causing extreme values.

We modeled the parameters $\mu,\sigma,k$ as suggested by [32] – assume
$\mu(x)=a+bx$, $\sigma(x)=e^{c+dx}$, $k(x)=e^{f+gx}$, where
$\\{a,b,c,d,f,g\\}$ is the new set of parameters. For each chosen parameter
set $\\{a,b,c,d,f,g\\}$, we formed $N$ bivariate samples of size $n$ by

  1. 1.

Sample $x_{1},x_{2},\dots,x_{n}$ independently from a uniform distribution;

  2. 2.

for all $i\in\\{1,2,\dots,n\\},$ generate $y_{i}$ from
$\operatorname*{GG}(\mu(x_{i}),\sigma(x_{i}),k(x_{i}))$ using one of the two
methods described in the previous paragraph.

We then applied to each sample MM/Double Kernel with a variety of covariate
knots and quantile order bases/bandwidths choices. Integrated Mean Squared
Error (IMSE), given in Equation 3.20, was calculated using the theoretical
quantile function Equation 4.10 at $q=0.1,0.25,0.5,0.75,0.9$ for each
approach. Estimated quantile curves based on the mean conditional predictions
$\hat{y}(x|q)=\hat{Q}_{q}(y|x)$ were also made to compare with the theoretical
ones.

We tried out four sets of parameter choices for $\\{a,b,c,d,f,g\\}$

  1. 1.

$a=5.39,b=-4.26,c=-2.52,d=3.25,f=-4.53,g=12.96$

  2. 2.

$a=1.384,b=0.092,c=-1.021,d=0.008,f=-3.493,g=4.766$

  3. 3.

$a=-1.5,b=-6,c=0.5,d=1.5,f=1.386,g=0$

  4. 4.

$a=-2,b=-0.75,c=-0.5,d=-4,f=-0.2,g=-1$,

among which the first two are approximate parameter estimates obtained by [32]
when fitting generalized gamma distributions to the aforementioned “WTD” data
set and an “IgG” data set [24] measuring concentrations of serum
immunoglobulin G for 298 kids under six, and the last two were argued by [32]
to produce particular interesting GG distribution scenarios. Theoretical
quantile curves at $q=0.1,0.25,0.5,0.75,0.9$ for the four generalized gamma
distributions are displayed in Figure 14.

IMSE charts produced using both Items 1 and 2 (Tables 6 and 7) largely favor
MM over Double Kernel (DK), as performances of DK with optimal bandwidth
choices are worse than those of most MM scenarios we have carried out. MM’s
superiority is also quite evidently showcased in the estimated quantile curves
at $q=0.1,0.25,0.5,0.75,0.9$ (Figures 15, 16 and 17 for Item 1 and Figures 18,
19 and 20 for Item 2), as MM’s plots not only resemble the theoretical ones
more but also smooth out rq’s empirical estimates (Figures B.20 and B.21) more
efficiently. We got similar results when using Item 4, with IMSE chart and
quantile curves presented in Table B.1 and Figures B.22, B.23, B.24 and B.25.
For Item 3, although the IMSE chart (Table B.2) leans towards DK with larger
$h_{1}$, double kernel frequently encounters extreme values when calculating
Equation 4.5 $\hat{Q}_{q}(y|x)=\hat{\hat{F}}^{-1}(q|x)$ and is thus unable to
produce estimated quantile curves. This is probably due to the manifest upward
shift of the theoretical quantile curve at $q=0.9$ (Figure 14(c)) for Item 3,
which significantly deviates from other quantile curves. MM, on the other
hand, does not run into such troubles. MM Estimated quantile curves
corresponding to Item 3 are presented in Figures B.26, B.27 and B.28.

(a) GG parameters as in Item 1

(b) GG parameters as in Item 2

(c) GG parameters as in Item 3

(d) GG parameters as in Item 4

Figure 14: Theoretical quantile curves at $q=0.1,0.25,0.5,0.75,0.9$ for
generalized gamma distributions with parameters specified in Items 1, 2, 3 and
4.

Quantile Level | 0.1 | 0.25 | 0.5 | 0.75 | 0.9
---|---|---|---|---|---
MM xseq3 | 15731.13211 | 18804.55677 | 20335.04031 | 19704.53243 | 13658.66566
MM xseq4 | 10826.00723 | 12304.97568 | 11508.04138 | 8792.64903 | 21247.54224
MM xseq5 | 10401.75120 | 10459.24631 | 9369.14093 | 8429.76193 | 22886.49413
MM xseq6 | 10966.40367 | 10287.35479 | 8902.93617 | 8590.76363 | 26394.91189
MM xseq7 | 11974.65885 | 10723.13314 | 9220.05222 | 9700.76951 | 27546.85202
MM xasym7 | 9474.81927 | 8978.33847 | 7269.27020 | 9559.08669 | 39670.68705
MM xseq3qseq3 | 16387.32192 | 15727.46428 | 20910.89332 | 23804.88056 | 13478.46414
MM xseq3qseq5 | 14561.38682 | 18765.77083 | 21887.80983 | 19337.62370 | 15196.28581
MM xseq4qseq3 | 11224.40731 | 9499.39642 | 12551.57089 | 8138.40265 | 25577.77186
MM xseq4qseq5 | 9859.05381 | 13058.06583 | 13016.67108 | 8032.73660 | 21423.01303
MM xseq5qseq3 | 11196.21847 | 8099.60148 | 10450.33535 | 7358.03035 | 28013.99754
MM xseq5qseq5 | 9335.08241 | 11684.91635 | 11032.72282 | 7620.14431 | 24441.25870
MM xseq6qseq3 | 12140.48855 | 8116.84493 | 10060.05862 | 7232.36086 | 31566.78552
MM xseq6qseq5 | 10130.85695 | 12129.81512 | 10594.15082 | 7924.11871 | 29482.06018
DK $h_{1}$=0.01 | 34904.68510 | 33209.31777 | 26198.68453 | 26869.14628 | 109294.90957
DK $h_{1}$=0.05 | 7090.81564 | 7420.21317 | 9356.54227 | 15685.48675 | 39204.01105
DK $h_{1}$=0.1 | 9375.15637 | 15014.34721 | 32748.95686 | 43978.01450 | 86078.69653
DK $h_{1}$=0.15 | 24540.11043 | 44337.87854 | 78907.50223 | 90945.24319 | 167363.10181
DK $h_{1}$=0.2 | 51385.92847 | 88727.01587 | 137268.44430 | 149175.56013 | 262179.51905
DK $h_{1}$=0.25 | 85067.72791 | 140664.69492 | 203164.43255 | 216404.63731 | 366280.14826
DK $h_{1}$=0.3 | 119018.90107 | 194034.30192 | 271070.24918 | 291794.60195 | 478713.31336

Table 6: IMSE (Equation 3.20) for MM and Double Kernel (DK) calculated using
$N=340$ Generalized Gamma samples of size $n=500$ with parameters specified in
Item 1 and $X\sim\operatorname*{Unif}(0,1)$. For MM, the first 6 settings
adopted logistic basis and the last 8 used natural spline basis (for quantile
orders) with equally spaced knots (numbers of knots for $q$ are specified
after “qseq”); we used equally spaced knots for the first 5 settings (numbers
of knots for $x$ are specified after “xseq”) and knots
$(0.2,0.4,0.5,0.55,0.6,0.7,0.9)$ for “MM xasym7” to natural-spline-transform
$x$. For DK, $h_{2}=10^{-4}$ for all choices of $h_{1}$.

For DK, $h_{1}=0.05$ is the clear best choice and $h_{1}=0.1$ the second best.
For MM, 4 and 5 knots for $x$ seem to be optimal, while 3 knots or asymmetric
knots for $x$ does not look adequate; the impact of quantile level basis
choices is likely more subtle. Comparing the two methods, DK with the best
bandwidth choice only outperforms MM with the worst choices, i.e., when we use
3 equally spaced knots or knots at $0.2,0.4,0.5,0.55,0.6,0.7,0.9$ to spline-
transform $x$.

Quantile Level | 0.1 | 0.25 | 0.5 | 0.75 | 0.9
---|---|---|---|---|---
MM xseq3 | 30.56947 | 33.33552 | 41.86750 | 52.84247 | 125.26318
MM xseq4 | 31.76553 | 31.64967 | 42.40184 | 60.49588 | 155.06402
MM xseq5 | 36.17384 | 32.25544 | 44.19438 | 69.03885 | 181.13795
MM xseq6 | 40.75397 | 35.59166 | 48.86846 | 78.87943 | 207.02464
MM xseq7 | 45.25154 | 39.47426 | 53.65562 | 86.63728 | 231.97000
MM xseq3qseq3 | 33.28659 | 32.54804 | 45.00762 | 89.05037 | 149.66787
MM xseq3qseq5 | 29.29923 | 35.40178 | 46.67796 | 65.92967 | 133.88833
MM xseq4qseq3 | 34.18627 | 29.94292 | 46.23849 | 98.24512 | 172.83624
MM xseq4qseq5 | 30.35540 | 33.10632 | 48.40875 | 76.73518 | 165.75680
MM xseq5qseq3 | 38.69695 | 30.26309 | 48.44576 | 108.71736 | 195.62337
MM xseq5qseq5 | 35.13168 | 33.81109 | 51.19494 | 87.02527 | 192.82714
MM xseq6qseq3 | 43.91158 | 33.24255 | 53.87782 | 121.20666 | 222.58875
MM xseq6qseq5 | 39.39483 | 37.49498 | 57.37674 | 100.79159 | 222.40250
DK $h_{1}$=0.06 | 133.45490 | 138.95991 | 196.71160 | 358.27465 | 888.78713
DK $h_{1}$=0.3 | 33.99643 | 39.25174 | 55.45666 | 83.58204 | 188.64835
DK $h_{1}$=0.6 | 42.30075 | 54.07889 | 73.11598 | 78.06503 | 146.32382
DK $h_{1}$=0.9 | 60.38152 | 89.56935 | 108.77301 | 107.50316 | 190.23795
DK $h_{1}$=1.2 | 81.53085 | 130.20296 | 148.07341 | 152.81118 | 273.73596
DK $h_{1}$=1.5 | 109.60298 | 172.80105 | 193.71430 | 214.61463 | 390.77361
DK $h_{1}$=1.8 | 144.71398 | 217.76392 | 245.52269 | 289.82433 | 535.43590

Table 7: IMSE (Equation 3.20) for MM and Double Kernel (DK) calculated using
$N=350$ Generalized Gamma samples of size $n=500$ with parameters specified in
Item 2 and $X\sim\operatorname*{Unif}(0,6)$. For MM, the first 5 settings
adopted logistic basis and the last 8 used natural spline basis (for quantile
orders) with equally spaced knots (numbers of knots for $q$ are specified
after “qseq”); we used equally spaced knots to natural-spline-transform $x$
throughout (numbers of knots for $x$ are specified after “xseq”). For DK,
$h_{2}=10^{-4}$ for all choices of $h_{1}$. Note that our choices for $h_{1}$
are larger since $X\sim\operatorname*{Unif}(0,6)$, not
$\operatorname*{Unif}(0,1)$ here.

For DK, the best two choices are still the second and third smallest –
$h_{1}=0.3$ and $h_{1}=0.6$, whose average performances are much closer now.
For MM, it seems this time that the fewer $x$ knots there are, the better the
algorithm; the impact of quantile-level basis choices is still more subtle,
but it seems (for this parameter set) that for a given $x$ transformation, the
logistic basis is the best of the three, and natural spline basis with 5 knots
outperforms that with 3 knots. MM still clearly beats DK on average, as DK
with optimal bandwidth choices only performs as well as MM with 5 knots for
$x$ (and significantly worse than MM with 3 or 4 $x$ knots).

Figure 15: Mean quantile curves estimated by Double Kernel for simulated GG
data with parameters as in Item 1 (theoretical curves – Figure 14(a)).
$h_{2}=10^{-4}$ for all cases, and $h_{1}=0.05,0.1,0.15,0.2,0.25,0.3$ from top
to bottom and left to right (the top two plots have $h_{1}=0.05,0.1$).

Figure 16: Mean quantile curves estimated by MM with logistic quantile level
basis for simulated GG data with parameters as in Item 1 (theoretical curves –
Figure 14(a)). The knots used to natural-spline-transform $x$ are $3,4,5,6,7$
equally spaced knots and $(0.2,0.4,0.5,0.55,0.6,0.7,0.9)$ from top to bottom
and left to right.

Figure 17: Mean quantile curves estimated by MM with natural spline (ns)
quantile level basis for simulated GG data with parameters as in Item 1
(theoretical curves – Figure 14(a)). The knots used to natural-spline-
transform $x$ are $3,4,5,6$ equally spaced knots from top to bottom. For each
row, the left/right plot used ns quantile level basis with 3/5 equally spaced
knots respectively.

Figure 18: Mean quantile curves estimated by Double Kernel for simulated GG
data with parameters as in Item 2 (theoretical curves – Figure 14(b)).
$h_{2}=10^{-4}$ for all cases, and $h_{1}=0.3,0.6,0.9,1.2,1.5,1.8$ from top to
bottom and left to right (the top two plots have $h_{1}=0.3,0.6$).

Figure 19: Mean quantile curves estimated by MM with logistic quantile level
basis for simulated GG data with parameters as in Item 2 (theoretical curves –
Figure 14(b)). The knots used to natural-spline-transform $x$ are $3,4,5,6,7$
equally spaced knots from top to bottom and left to right.

Figure 20: Mean quantile curves estimated by MM with natural spline (ns)
quantile level basis for simulated GG data with parameters as in Item 2
(theoretical curves – Figure 14(b)). The knots used to natural-spline-
transform $x$ are $3,4,5,6$ equally spaced knots from top to bottom. For each
row, the left/right plot used ns quantile level basis with 3/5 equally spaced
knots, respectively.

##  5 Regularized and Monotone Quantile Regression

We now discuss practical techniques commonly adopted in two quantile
regression settings with special requirements – adaptive lasso penalty for
regularized regression and monotone B-splines for monotone regression.

###  5.1 Adaptive Lasso for Heterogeneous Variable Selection

In diverse sparse modeling settings requiring efficient predictor selection,
the adaptive lasso has been proven highly competent in consistent variable
selection by imposing an adequately weighted $\mathit{l}_{1}$ penalty
$\lambda\sum_{j=1}^{p}w_{j}|\beta_{j}|,\lambda\geq 0$ [48]. Several authors
(e.g., [45], [33]) have extended this powerful technique to account for
possible quantile-varying covariate effects by defining a penalized quantile
regression problem with parameter estimate202020We write
$\bm{\beta}=(\beta_{1},\beta_{2},\dots,\beta_{p})^{T}$, where $\beta_{1}$
denotes the coefficient for the intercept term.

$\displaystyle\hat{\bm{\beta}}_{\lambda}(q)_{p\times 1}=$
$\displaystyle\underset{\bm{\beta}}{\text{ argmin}}\;L_{q}(\bm{\beta}),\text{
where}$ $\displaystyle L_{q}(\bm{\beta})=$
$\displaystyle\sum_{i=1}^{n}\rho_{q}(r_{i}(\bm{\beta}))+\lambda\sum_{j=2}^{p}\frac{|\beta_{j}|}{\big{|}\tilde{\beta}_{j}(q)\big{|}}\text{
and}$ (5.1) $\displaystyle\tilde{\bm{\beta}}(q)_{p\times 1}=$
$\displaystyle\underset{\bm{\beta}}{\text{
argmin}}\;\sum_{i=1}^{n}\rho_{q}(r_{i}(\bm{\beta}))\text{ (the nonregularized
quantile regression estimator)}$

at any quantile level $q\in(0,1)$.

Thanks to the AM-GM Inequality,

$Q_{q,\operatorname*{lasso}}\left(\bm{\beta}\bigg{|}\bm{\beta}^{(t)}\right)=\sum_{i=1}^{n}\zeta_{q}\left(r_{i}\bigg{|}r_{i}^{(t)}\right)+\lambda\sum_{j=2}^{p}\frac{1}{\left|\tilde{\beta}_{j}(q)\right|}\left[\frac{1}{2}\left|\beta_{j}^{(t)}\right|+\frac{\beta_{j}^{2}}{2\left|\beta_{j}^{(t)}\right|}\right]$
(5.2)

is a majorizer of Section 5.1 at $\bm{\beta}^{(t)}$. To mend the flaw that
Equation 5.2 is undefined when $r_{i}^{(t)}=0$ for any
$i\in\\{1,2,\ldots,n\\}$ or $\beta_{j}^{(t)}=0$ for any
$j\in\\{2,3,\dots,p\\}$, a realistic concern as
$r_{i}\big{(}\bm{\beta}^{(t)}\big{)}$, $i\in\\{1,2,\ldots,n\\}$ for some
$f_{i}\big{(}\bm{\beta}^{(t)}\big{)}$s and $\beta_{j}^{(t)}$,
$j\in\\{2,3,\dots,p\\}$ for penalized quantile regression may well converge to
0 as $t\rightarrow\infty$, we can adopt212121Note that while it is still
ensured that $L_{q}^{\epsilon}(\bm{\beta})\leq
Q^{\epsilon_{l}}_{q,\operatorname*{lasso}}(\bm{\beta}|\bm{\beta}^{(t)})\;\forall\;\bm{\beta}$,
introducing a perturbation $\epsilon_{l}>0$ makes
$L_{q}^{\epsilon}(\bm{\beta}^{(t)})\neq
Q^{\epsilon_{l}}_{q,\operatorname*{lasso}}(\bm{\beta}^{(t)}|\bm{\beta}^{(t)})$.

$\displaystyle
Q^{\epsilon_{l}}_{q,\operatorname*{lasso}}\left(\bm{\beta}\bigg{|}\bm{\beta}^{(t)}\right)=\sum_{i=1}^{n}\zeta_{q}^{\epsilon}\left(r_{i}\bigg{|}r_{i}^{(t)}\right)+\lambda\sum_{j=2}^{p}\frac{1}{\big{|}\tilde{\beta}_{j}(q)\big{|}}\left[\frac{1}{2}\left(\big{|}\beta_{j}^{(t)}\big{|}+\epsilon_{l}\right)+\frac{\beta_{j}^{2}}{2\left(\big{|}\beta_{j}^{(t)}\big{|}+\epsilon_{l}\right)}\right],$
(5.3)

a majorizer at $\bm{\beta}^{(t)}$ of

$L_{q}^{\epsilon}(\bm{\beta})=\sum_{i=1}^{n}\rho_{q}^{\epsilon}(r_{i}(\bm{\beta}))+\lambda\sum_{j=2}^{p}\frac{|\beta_{j}|}{|\tilde{\beta}_{j}(q)|}=\sum_{i=1}^{n}\left\\{\rho_{q}(r_{i}(\bm{\beta}))-\frac{\epsilon}{2}\ln(\epsilon+|r_{i}(\bm{\beta})|)\right\\}+\lambda\sum_{j=2}^{p}\frac{|\beta_{j}|}{|\tilde{\beta}_{j}(q)|},$
(5.4)

which closely approximates $L_{q}(\bm{\beta})$ in Section 5.1. Let

$\displaystyle\bm{V}^{(t)}_{p\times p}$ $\displaystyle=\begin{pmatrix}0&&&\\\
&\frac{1}{|\tilde{\beta}_{2}(q)|\cdot(|\beta_{2}^{(t)}|+\epsilon_{l})}&&\\\
&&\ddots&\\\
&&&\frac{1}{|\tilde{\beta}_{p}(q)|\cdot(|\beta_{p}^{(t)}|+\epsilon_{l})}\end{pmatrix},$
(5.5) $\displaystyle\bm{W}^{(t)}_{n\times n}$
$\displaystyle=\begin{pmatrix}\frac{1}{\epsilon+|r_{1}(\bm{\beta}^{(t)})|}&&&\\\
&\frac{1}{\epsilon+|r_{2}(\bm{\beta}^{(t)})|}&&\\\ &&\ddots&\\\
&&&\frac{1}{\epsilon+|r_{n}(\bm{\beta}^{(t)})|}\end{pmatrix},\text{ and }$
$\displaystyle\bm{c}_{n\times 1}$ $\displaystyle=(4q-2,4q-2,\dots,4q-2)^{T}.$

Then

$Q^{\epsilon_{l}}_{q,\operatorname*{lasso}}(\bm{\beta}|\bm{\beta}^{(t)})=\frac{1}{4}\left(\bm{(y-X\beta)}^{T}\bm{W}^{(t)}\bm{(y-X\beta)}+\bm{c}^{T}(\bm{y-X\beta})\right)+\frac{\lambda}{2}\bm{\beta}^{T}\bm{V}^{(t)}\bm{\beta}+\operatorname*{const}.$
(5.6)

Hence, minimizing
$Q^{\epsilon_{l}}_{q,\operatorname*{lasso}}(\bm{\beta}|\bm{\beta}^{(t)})$ is
equivalent to minimizing

$\displaystyle h(\bm{\beta}|\bm{\beta}^{(t)})$
$\displaystyle=g(\bm{\beta}|\bm{\beta}^{(t)})+2\lambda\bm{\beta}^{T}\bm{V}^{(t)}\bm{\beta},\text{
where }$ (5.7) $\displaystyle g(\bm{\beta}|\bm{\beta}^{(t)})$
$\displaystyle=\bm{(y-X\beta)}^{T}\bm{W}^{(t)}\bm{(y-X\beta)}+\bm{c}^{T}(\bm{y-X\beta}).$

We have shown in Section 3.1 that $\bm{X}^{T}\bm{W}^{(t)}\bm{X}$ is positive
definite and thus invertible for a full rank $\bm{X}_{n\times p}$. It is also
clear that $\bm{V}^{(t)}$ is positive semi-definite, as it is a diagonal
matrix whose entries are all non-negative. Hence, for a full rank
$\bm{X}_{n\times p}$, the Hessian

$\displaystyle\frac{\partial^{2}}{\partial\bm{\beta}\partial\bm{\beta}^{T}}h(\bm{\beta}|\bm{\beta}^{(t)})=2\bm{X}^{T}\bm{W}^{(t)}\bm{X}+4\lambda\bm{V}^{(t)}\;\text{
is positive definite since }\lambda\geq 0.$ (5.8)

We can thus proceed to solve for $\bm{\beta}^{(t+1)}$ by directly minimizing
Equation 5.7:

$\displaystyle\frac{\partial}{\partial\bm{\beta}}h(\bm{\beta}|\bm{\beta}^{(t)})\bigg{\rvert}_{\bm{\beta}=\bm{\beta}^{(t+1)}}$
$\displaystyle=\left(2\bm{X}^{T}\bm{W}^{(t)}\bm{X\beta}-2\bm{X}^{T}\bm{W}^{(t)}\bm{y}-\bm{X}^{T}\bm{c}+4\lambda\bm{V}^{(t)}\bm{\beta}\right)\bigg{\rvert}_{\bm{\beta}=\bm{\beta}^{(t+1)}}$
$\displaystyle=2\left(\bm{X}^{T}\bm{W}^{(t)}\bm{X}+2\lambda\bm{V}^{(t)}\right)\bm{\beta}^{(t+1)}-2\bm{X}^{T}\bm{W}^{(t)}\bm{y}-\bm{X}^{T}\bm{c}=\bm{0}$
$\displaystyle\Rightarrow\bm{\beta}^{(t+1)}=\frac{1}{2}$
$\displaystyle\left(\bm{X}^{T}\bm{W}^{(t)}\bm{X}+2\lambda\bm{V}^{(t)}\right)^{-1}\bm{X}^{T}\left[2\bm{W}^{(t)}\bm{y}+\bm{c}\right].$
(5.9)

$\left(\bm{X}^{T}\bm{W}^{(t)}\bm{X}+2\lambda\bm{V}^{(t)}\right)^{-1}$ indeed
exists for a full rank $\bm{X}$ as Equation 5.8 $\Rightarrow
2\bm{X}^{T}\bm{W}^{(t)}\bm{X}+4\lambda\bm{V}^{(t)}=2(\bm{X}^{T}\bm{W}^{(t)}\bm{X}+2\lambda\bm{V}^{(t)})\text{
is invertible }\iff\bm{X}^{T}\bm{W}^{(t)}\bm{X}+2\lambda\bm{V}^{(t)}$ is
invertible.

It remains to find suitable choices of $\lambda>0$ for each $q\in(0,1)$ of
interest. We followed [33]’s idea and used a quantile regression BIC metric
put forward by [30]:

${\operatorname*{BIC}}_{q}(\lambda)={\operatorname*{BIC}}_{q}\left(\hat{\bm{\beta}}_{\lambda}(q)\right)=\ln\left(\sum_{i=1}^{n}\rho_{q}\left(r_{i}\left(\hat{\bm{\beta}}_{\lambda}(q)\right)\right)\right)+|\mathit{S}_{q}^{\lambda}|\frac{\ln(n)}{2n},$
(5.10)

where $|\mathit{S}_{q}^{\lambda}|$ represents the size of
$\mathit{S}_{q}^{\lambda}=\\{2\leq j\leq
p:|(\hat{\beta}_{\lambda})_{j}|>0\\}$.222222When actually implementing the
algorithm in statistical software like R, we set a threshold and take
$(\hat{\beta}_{\lambda})_{j}=0$ if $|(\hat{\beta}_{\lambda})_{j}|$ is less
than this threshold value. A derivation to Equation 5.10 is presented in
Appendix A Proposition A.2. Our algorithm thus proceeds according to the
following steps:

  1. 1.

Choose a quantile level $q\in(0,1)$ of interest;

  2. 2.

Calculate a nonregularized estimator $\tilde{\bm{\beta}}(q)$ (using rq or MM);

  3. 3.

Propose a sequence of candidate $\lambda\in(0,1)$ choices;

  4. 4.

For each candidate $\lambda$, obtain a $\hat{\bm{\beta}}_{\lambda}(q)$ that
attempts to minimize $L_{q}^{\epsilon}(\bm{\beta})$ in Equation 5.4 via our MM
algorithm detailed above;

  5. 5.

Evaluate ${\operatorname*{BIC}}_{q}(\lambda)$ in Equation 5.10 for this
coefficient estimate $\hat{\bm{\beta}}_{\lambda}(q)$;

  6. 6.

Find $\lambda_{\operatorname*{opt}}(q)$ among our picked sequence of candidate
$\lambda$’s that yields the smallest ${\operatorname*{BIC}}_{q}(\lambda)$;

  7. 7.

Calculate our regularized estimator
$\hat{\bm{\beta}}_{\lambda_{\operatorname*{opt}}}(q)$ with this
$\lambda_{\operatorname*{opt}}$ at the desired quantile level $q$.

$\bm{q=0.25}$ | $\bm{q=0.5}$ | $\bm{q=0.75}$
---|---|---
$x_{j}$ | $(\hat{\beta}_{\lambda_{\operatorname*{opt}}})_{j}$ | $x_{j}$ | $(\hat{\beta}_{\lambda_{\operatorname*{opt}}})_{j}$ | $x_{j}$ | $(\hat{\beta}_{\lambda_{\operatorname*{opt}}})_{j}$
$x_{9}$ | 0.05258 | $x_{12}$ | -0.03797 | $x_{12}$ | -0.04276
$x_{14}$ | 0.02368 | $x_{9}$ | 0.03291 | $x_{13}$ | 0.03740
$x_{3}$ | -0.02063 | $x_{13}$ | 0.03278 | $x_{9}$ | 0.03614
$x_{10}$ | -0.01457 | $x_{1}$ | 0.01673 | $x_{6}$ | -0.02074
$x_{7}$ | -0.01440 | $x_{14}$ | 0.01394 | $x_{2}$ | -0.01655
$x_{5}$ | -0.01329 | $x_{8}$ | 0.01008 | $x_{8}$ | 0.00873
$x_{1}$ | 0.01167 | $x_{10}$ | -0.00716 | $x_{4}$ | 5.76206E-10
$x_{11}$ | -0.01087 | $x_{2}$ | -0.00308 | $x_{14}$ | 4.89131E-10
$x_{2}$ | -0.00938 | $x_{3}$ | -0.00300 | $x_{7}$ | -4.87400E-10
$x_{13}$ | -0.00933 | $x_{5}$ | -1.07803E-09 | $x_{1}$ | 4.21575E-10
$x_{15}$ | -0.00780 | $x_{4}$ | 3.95779E-10 | $x_{5}$ | -4.98012E-11
$x_{12}$ | -1.09786E-10 | $x_{11}$ | -3.90428E-11 | $x_{15}$ | 4.23412E-11
$x_{8}$ | -1.58700E-12 | $x_{15}$ | 1.70300E-11 | $x_{10}$ | 3.32322E-11
$x_{4}$ | 1.53641E-12 | $x_{6}$ | -1.29810E-11 | $x_{3}$ | -2.90341E-12
$x_{6}$ | 1.36344E-13 | $x_{7}$ | -9.11257E-12 | $x_{11}$ | -1.24519E-12
Table 8: Coefficient estimates obtained via MM in Section 5.1 for regularized
quantile regression with adaptive lasso penalty at $q=0.25,0.5,0.75$ on the
“pollution” data set. For each $q$, the covariates and their respective
penalized coefficient estimates are presented in decreasing order of the
magnitude of the absolute values of the estimators.

We started off with the “pollution” data set [31] first raised in Section 3.1,
as [33] have produced some “guideline” results obtained by their adaptive
lasso “local” quantile regression methods $SS(q)$ at $q=0.25,0.5,0.75$. As
presented in Table 8, our results obtained from MM do align quite closely with
[33]’s counterparts (Table 9). We gained more confidence in our MM algorithm
by further comparing our MM penalized estimates with those from quantreg’s
rq.fit.lasso (Table 10) with the penalty vector specified as
$\left(0,\frac{\hat{\lambda}_{\operatorname*{opt}}(q)}{|\tilde{\beta}_{2}(q)|},\dots,\frac{\hat{\lambda}_{\operatorname*{opt}}(q)}{|\tilde{\beta}_{p}(q)|}\right)$,
where $\hat{\lambda}_{\operatorname*{opt}}(q)$ is the $\lambda\in(0,1)$ we
found that yields the smallest $\operatorname*{BIC}_{q}$.

Method | Selected Variables
---|---
$SS(0.25)$ | $x_{1},x_{3},x_{9},x_{10},x_{13},x_{14}$
$SS(0.5)$ | $x_{1},x_{8},x_{9},x_{10},x_{12},x_{14}$
$SS(0.75)$ | $x_{6},x_{8},x_{9},x_{12},x_{14}$
Table 9: Variables (of the “pollution” data set) selected by penalized
quantile regression methods with adaptive lasso penalty [33].

Effect | MM | rq.fit.lasso | MM | rq.fit.lasso | MM | rq.fit.lasso
---|---|---|---|---|---|---
(intercept) | 6.82237 | 6.82237 | 6.84400 | 6.84499 | 6.87363 | 6.86302
prec | 0.01167 | 0.01167 | 0.01673 | 0.01735 | 4.22E-10 | 0.01056
jant | -0.00938 | -0.00938 | -0.00308 | -0.00557 | -0.01655 | -0.01672
jult | -0.02063 | -0.02063 | -0.00300 | -0.00237 | -2.90E-12 | -8.74E-11
ovr95 | 1.54E-12 | 7.20E-15 | 3.96E-10 | 1.23E-12 | 5.76E-10 | 3.32E-09
popn | -0.01329 | -0.01329 | -1.08E-09 | -0.00153 | -4.98E-11 | -0.00389
educ | 1.36E-13 | 1.34E-15 | -1.30E-11 | -3.51E-14 | -0.02074 | -0.01715
hous | -0.01440 | -0.01440 | -9.11E-12 | -2.50E-14 | -4.87E-10 | -0.00362
dens | -1.59E-12 | -1.03E-14 | 0.01008 | 0.00948 | 0.00873 | 0.00883
nonw | 0.05258 | 0.05258 | 0.03291 | 0.03455 | 0.03614 | 0.04168
wwdrk | -0.01457 | -0.01457 | -0.00716 | -0.00814 | 3.32E-11 | 4.11E-10
poor | -0.01087 | -0.01087 | -3.90E-11 | 7.51E-13 | -1.25E-12 | -0.00177
hc | -1.10E-10 | -1.63E-12 | -0.03797 | -0.03740 | -0.04276 | -0.05761
nox | -0.00933 | -0.00933 | 0.03278 | 0.03298 | 0.03740 | 0.05906
so | 0.02368 | 0.02368 | 0.01394 | 0.01366 | 4.89E-10 | 7.30E-10
humid | -0.00780 | -0.00780 | 1.70E-11 | -1.09E-14 | 4.23E-11 | 0.00272

Table 10: Lasso penalized quantile regression coefficient estimates obtained
from MM and rq.fit.lasso for the “pollution” data set at $q=0.25,0.5,0.75$.

We then implemented our MM algorithm on three other high-dimensional data sets
at $q=0.1,0.25,0.5,0.75,0.9$ to uncover potential non-constant predictor
effects in diverse fields.

Economics Data

We constructed our first data set by compiling US monthly time series data for
9 economic variables (unemployment rate $y$, personal consumption expenditures
(PCE) $x_{1}$, personal savings rate $x_{2}$, population $x_{3}$, Industrial
Production Total Index $x_{4}$, Consumer Price Index (CPI) $x_{5}$, Policy
Uncertainty Index232323Based on US News. $x_{6}$, Producer Price Index by
Commodity $x_{9}$, and KC Fed Labor Market Conditions Index $x_{10}$)
downloaded from [13] and US monthly Trade Policy Uncertainty (TPU),
Geopolitical Risk (GPR) Indices $x_{7},x_{8}$ constructed by [5] and [4]. We
considered observations (historical months) that have data in all 11 variables
depicted above, and hence our finalized “econ” data set comprises $n=347$
continuous months from January 1992 to November 2020 and $p-1=10$ explanatory
variables.

$\bm{q=0.1}$ | $\bm{q=0.25}$ | $\bm{q=0.5}$ | $\bm{q=0.75}$ | $\bm{q=0.9}$
---|---|---|---|---
$x_{j}$ | $(\hat{\beta}_{\lambda_{\operatorname*{opt}}})_{j}$ | $x_{j}$ | $(\hat{\beta}_{\lambda_{\operatorname*{opt}}})_{j}$ | $x_{j}$ | $(\hat{\beta}_{\lambda_{\operatorname*{opt}}})_{j}$ | $x_{j}$ | $(\hat{\beta}_{\lambda_{\operatorname*{opt}}})_{j}$ | $x_{j}$ | $(\hat{\beta}_{\lambda_{\operatorname*{opt}}})_{j}$
$x_{1}$ | -6.07027 | $x_{3}$ | 9.06328 | $x_{3}$ | 11.07810 | $x_{3}$ | 10.27440 | $x_{3}$ | 8.19060
$x_{3}$ | 5.92319 | $x_{1}$ | -5.26884 | $x_{5}$ | -10.24142 | $x_{5}$ | -6.53847 | $x_{5}$ | -3.20889
$x_{9}$ | 1.62627 | $x_{5}$ | -4.93119 | $x_{9}$ | 3.14862 | $x_{4}$ | -3.07802 | $x_{4}$ | -3.15706
$x_{4}$ | -1.57856 | $x_{9}$ | 2.83706 | $x_{4}$ | -2.62308 | $x_{1}$ | -2.40785 | $x_{1}$ | -2.73452
$x_{8}$ | 0.00010 | $x_{4}$ | -1.89304 | $x_{1}$ | -1.43221 | $x_{9}$ | 1.77553 | $x_{9}$ | 0.97606
$x_{7}$ | 2.53E-09 | $x_{6}$ | 0.17831 | $x_{6}$ | 0.36607 | $x_{7}$ | -0.36456 | $x_{7}$ | -0.44530
$x_{10}$ | 1.55E-09 | $x_{8}$ | 0.09129 | $x_{10}$ | -0.24760 | $x_{6}$ | 0.16250 | $x_{6}$ | 0.21889
$x_{2}$ | 9.98E-10 | $x_{7}$ | 7.52E-07 | $x_{7}$ | -0.19279 | $x_{2}$ | 0.15769 | $x_{2}$ | 0.14194
$x_{5}$ | -6.33E-10 | $x_{10}$ | -6.26E-09 | $x_{8}$ | 1.02E-06 | $x_{10}$ | -0.14640 | $x_{10}$ | -2.99E-07
$x_{6}$ | 1.00E-11 | $x_{2}$ | -1.70E-11 | $x_{2}$ | 1.48E-08 | $x_{8}$ | 3.72E-10 | $x_{8}$ | -8.40E-08
Table 11: Coefficient estimates obtained via MM in Section 5.1 for regularized
quantile regression with adaptive lasso penalty at $q=0.1,0.25,0.5,0.75,0.9$
on the “econ” data set. For each $q$, the covariates and their respective
penalized coefficient estimates are presented in decreasing order of the
magnitude of the absolute values of the estimators.

Our results, as shown in Table 11, lend credence to non-constant predictor
effects. Although the highest correlated predictor population $x_{3}$’s
effects (positively influencing $y$) are relatively constant, 2 other critical
predictors showcase manifest non-constant effects across various quantile
orders: CPI $x_{5}$ is insignificant at $q=0.1$ yet counts a lot at higher
quantiles; PCE $x_{1}$ appears clearly more significant (in negatively
affecting $y$) at lower quantiles. Our finding that the inflation rate
(measured by CPI $x_{5}$) negatively correlates to unemployment rate $y$
significantly when $y$ is larger aligns closely to the current widely accepted
expectations-augmented Philips curve theory.242424That cyclical unemployment
and unanticipated inflation are inversely related. In addition, $x_{4}$
(negative correlation with $y$) seems to matter more at higher $q$ and $x_{9}$
(positive correlation with $y$) at lower $q$. These insightful findings, if
backed by more evidence, may aid governments considerably in coming up with
more effective and adequate macroeconomic policies to combat high unemployment
rates in diverse situations.

Social Data

Our second data set is a “social” one obtained by combining certain columns of
a public data set investigating worldwide country-level life expectancy [43]
and other potential predictors’ per-country data gathered from [42], [41],
[44], and [40]. We removed variables with insufficient non-missing data and
considered $n=111$ countries with no missing data on all remaining 21 (1
explanatory variable $y$, Life Expectancy, and 20 explanatory
variables252525Hence, $p=20+1=21$ (including the intercept). The 20
coviariates are crude suicide rate $x_{1}$, Corruption Perceptions Index (CPI)
$x_{2}$, fertility rate $x_{3}$, mobilization resources $x_{4}$, unemployment
rate $x_{5}$, adult mortality rate $x_{6}$, infant death frequency $x_{7}$,
Hepatitis B immunization coverage $x_{8}$, measles frequency $x_{9}$, average
Body Mass Index (BMI) $x_{10}$, under-five deaths $x_{11}$, polio and DTP3
immunization coverages $x_{12}$ and $x_{13}$, HIV/AIDS deaths $x_{14}$, GDP
per capita $x_{15}$, population $x_{16}$, thinness frequency among 10-19-year-
olds and among 5-9-year-olds $x_{17}$ and $x_{18}$, income composition of
resources $x_{19}$, and number of school years $x_{20}$.) variables.

Should there be drastic non-constant predictor effects, in which case
regularized quantile regression may greatly improve fit, a country’s leaders
can then devise more potent policies to influence its life expectancy in its
desired way given the quantile of life expectancy the country stands
worldwide, by controlling the corresponding influential variables around that
$q$.

$\bm{q=0.1}$ | $\bm{q=0.25}$ | $\bm{q=0.5}$ | $\bm{q=0.75}$ | $\bm{q=0.9}$
---|---|---|---|---
$x_{j}$ | $(\hat{\beta}_{\lambda_{\operatorname*{opt}}})_{j}$ | $x_{j}$ | $(\hat{\beta}_{\lambda_{\operatorname*{opt}}})_{j}$ | $x_{j}$ | $(\hat{\beta}_{\lambda_{\operatorname*{opt}}})_{j}$ | $x_{j}$ | $(\hat{\beta}_{\lambda_{\operatorname*{opt}}})_{j}$ | $x_{j}$ | $(\hat{\beta}_{\lambda_{\operatorname*{opt}}})_{j}$
$x_{19}$ | 3.72627 | $x_{3}$ | -3.11104 | $x_{19}$ | 3.09364 | $x_{19}$ | 3.79070 | $x_{6}$ | -3.43217
$x_{6}$ | -2.44059 | $x_{19}$ | 2.23081 | $x_{6}$ | -2.54648 | $x_{6}$ | -2.00230 | $x_{19}$ | 3.20339
$x_{3}$ | -2.38420 | $x_{6}$ | -1.71024 | $x_{3}$ | -1.96883 | $x_{3}$ | -1.50224 | $x_{3}$ | -1.03290
$x_{8}$ | 1.36639 | $x_{1}$ | -0.81732 | $x_{1}$ | -0.87842 | $x_{1}$ | -0.96843 | $x_{14}$ | 0.98214
$x_{17}$ | -1.00701 | $x_{2}$ | 0.79934 | $x_{4}$ | 0.53638 | $x_{16}$ | -0.56308 | $x_{1}$ | -0.76843
$x_{11}$ | 0.85944 | $x_{17}$ | -0.77932 | $x_{16}$ | -0.43499 | $x_{8}$ | 0.54486 | $x_{8}$ | 0.57012
$x_{2}$ | 0.69624 | $x_{4}$ | 0.72324 | $x_{5}$ | -0.41553 | $x_{4}$ | 0.40565 | $x_{4}$ | 0.40426
$x_{13}$ | -0.63655 | $x_{8}$ | 0.67395 | $x_{2}$ | 0.30917 | $x_{2}$ | 0.26818 | $x_{18}$ | -0.34016
$x_{1}$ | -0.57268 | $x_{7}$ | 0.45584 | $x_{14}$ | -0.30400 | $x_{18}$ | -0.15954 | $x_{11}$ | -0.01845
$x_{14}$ | -0.50570 | $x_{16}$ | -0.31903 | $x_{8}$ | 0.15220 | $x_{5}$ | -0.13403 | $x_{16}$ | -1.88E-06
$x_{4}$ | 0.27249 | $x_{14}$ | -0.29717 | $x_{17}$ | -0.14157 | $x_{11}$ | -0.05771 | $x_{5}$ | 6.12E-07
$x_{16}$ | -0.17622 | $x_{5}$ | -0.29463 | $x_{12}$ | -0.00092 | $x_{12}$ | -2.12E-06 | $x_{10}$ | 2.26E-08
$x_{10}$ | -7.54E-07 | $x_{9}$ | -8.52E-05 | $x_{7}$ | 6.43E-08 | $x_{15}$ | -1.82E-08 | $x_{20}$ | -7.03E-09
$x_{18}$ | 6.10E-08 | $x_{11}$ | 1.90E-08 | $x_{9}$ | 2.46E-08 | $x_{7}$ | -3.43E-09 | $x_{7}$ | -4.48E-09
$x_{5}$ | -1.43E-08 | $x_{15}$ | 1.29E-08 | $x_{10}$ | 1.61E-08 | $x_{17}$ | 2.37E-09 | $x_{2}$ | 3.94E-09
$x_{9}$ | -5.34E-09 | $x_{13}$ | -3.60E-09 | $x_{18}$ | -7.17E-09 | $x_{10}$ | 1.76E-09 | $x_{15}$ | 2.58E-09
$x_{7}$ | -2.98E-09 | $x_{12}$ | 1.52E-09 | $x_{11}$ | 6.72E-09 | $x_{9}$ | -1.52E-09 | $x_{12}$ | 1.55E-09
$x_{12}$ | -1.38E-09 | $x_{18}$ | 2.20E-10 | $x_{20}$ | -6.27E-09 | $x_{20}$ | 6.87E-10 | $x_{9}$ | 1.16E-09
$x_{20}$ | 1.93E-10 | $x_{10}$ | 1.90E-10 | $x_{13}$ | -2.19E-09 | $x_{14}$ | -4.76E-11 | $x_{17}$ | -1.88E-10
$x_{15}$ | 2.80E-11 | $x_{20}$ | 3.49E-11 | $x_{15}$ | 7.01E-10 | $x_{13}$ | 3.86E-12 | $x_{13}$ | 1.48E-10
Table 12: Coefficient estimates obtained via MM in Section 5.1 for regularized
quantile regression with adaptive lasso penalty at $q=0.1,0.25,0.5,0.75,0.9$
on the “social” data set. For each $q$, the covariates and their respective
penalized coefficient estimates are presented in decreasing order of the
magnitude of the absolute values of the estimators.

It turns out that covariates for this “social” data set do not seem to exhibit
notably varying effects for different quantile orders (Table 12). The most
attributable predictors ($x_{19}$ – positively affects $y$; $x_{6},x_{1}$ –
negatively correlate to $y$) seem to possess quite constant effects across
quantiles. There are, however, covariates among second-tier influential
predictors that do suggest possible varying effects. Effects of $x_{17}$
(which affects $y$ negatively) and $x_{2}$ (which influences $y$ positively)
appear to decrease as $q$ increases, and $x_{8}$ (which positively correlates
to $y$) seems less important in middle quantiles.

Forest Fire Data

Our last data set is a real-world “forest fire” data set collected from
Northeast Portugal [10], which contains a response variable “burned area” as
well as potential spatial ($x_{1},x_{2}$: x, y coordinates on the Montesinho
park map), temporal ($x_{3}$: month of the year; $x_{4}$: day of the week),
meteorological weather ($x_{9}$: temperature; $x_{10}$: relative humidity;
$x_{11}$: wind speed; $x_{12}$: rainfall) explanatory variables and key Fire
Weather Index (FWI) components including the Fine Fuel Moisture Code (FFMC)
$x_{5}$, the Duff Moisture Code (DMC) $x_{6}$, the Drought Code (DC) $x_{7}$,
and the Initial Spread Index (ISI) $x_{8}$. We used a subset (270 of the 517
instances when there is a forest fire) of the original data set for analysis
and transformed the response variable using $y=\ln(x+0.1)$, as the original
response variable “area” is highly concentrated at $0_{+}$.

As small fires are less damaging but much more frequent and large fires are
less frequent but much more destructive, it would significantly aid targeted
firefighting and thus reduce fires’ negative impacts if we can discern
potential differentiated predictor effects across different fire scales and
make more accurate predictions on a given initiated fire’s possible damages.

$\bm{q=0.1}$ | $\bm{q=0.25}$ | $\bm{q=0.5}$ | $\bm{q=0.75}$ | $\bm{q=0.9}$
---|---|---|---|---
$x_{j}$ | $(\hat{\beta}_{\lambda_{\operatorname*{opt}}})_{j}$ | $x_{j}$ | $(\hat{\beta}_{\lambda_{\operatorname*{opt}}})_{j}$ | $x_{j}$ | $(\hat{\beta}_{\lambda_{\operatorname*{opt}}})_{j}$ | $x_{j}$ | $(\hat{\beta}_{\lambda_{\operatorname*{opt}}})_{j}$ | $x_{j}$ | $(\hat{\beta}_{\lambda_{\operatorname*{opt}}})_{j}$
$x_{9}$ | -0.22187 | $x_{9}$ | -0.18948 | $x_{7}$ | -0.15219 | $x_{7}$ | -0.58951 | $x_{9}$ | 0.23867
$x_{8}$ | -0.10141 | $x_{6}$ | 0.15076 | $x_{3}$ | 0.10327 | $x_{9}$ | 0.43749 | $x_{6}$ | 0.12303
$x_{4}$ | 0.05433 | $x_{5}$ | -0.11819 | $x_{10}$ | -0.03382 | $x_{3}$ | 0.37752 | $x_{1}$ | 0.07108
$x_{6}$ | 1.65E-07 | $x_{7}$ | -0.06580 | $x_{9}$ | -0.00555 | $x_{6}$ | 0.34871 | $x_{8}$ | -0.04887
$x_{12}$ | 7.54E-09 | $x_{1}$ | -0.00619 | $x_{6}$ | 3.11E-07 | $x_{2}$ | -0.32366 | $x_{7}$ | -2.11E-06
$x_{1}$ | 1.68E-09 | $x_{3}$ | 1.35E-07 | $x_{2}$ | -1.74E-08 | $x_{1}$ | 0.322171 | $x_{11}$ | 3.32E-08
$x_{10}$ | -1.19E-09 | $x_{10}$ | -1.70E-08 | $x_{1}$ | 1.60E-08 | $x_{8}$ | -0.26285 | $x_{3}$ | -8.89E-09
$x_{7}$ | -1.10E-09 | $x_{4}$ | 3.60E-09 | $x_{4}$ | 1.26E-08 | $x_{4}$ | 0.18100 | $x_{2}$ | -4.23E-09
$x_{3}$ | 5.74E-10 | $x_{12}$ | 3.13E-09 | $x_{12}$ | 7.57E-09 | $x_{11}$ | 0.14862 | $x_{12}$ | -3.95E-09
$x_{11}$ | 5.26E-10 | $x_{11}$ | 1.80E-09 | $x_{8}$ | -5.66E-09 | $x_{5}$ | 1.76E-06 | $x_{5}$ | -1.18E-09
$x_{5}$ | -2.56E-10 | $x_{8}$ | -1.19E-09 | $x_{5}$ | -1.75E-09 | $x_{12}$ | -3.12E-08 | $x_{10}$ | -9.88E-10
$x_{2}$ | 6.10E-11 | $x_{2}$ | 1.09E-09 | $x_{11}$ | 6.05E-10 | $x_{10}$ | -2.93E-09 | $x_{4}$ | 4.72E-10
Table 13: Coefficient estimates obtained via MM in Section 5.1 for regularized
quantile regression with adaptive lasso penalty at $q=0.1,0.25,0.5,0.75,0.9$
on the “forest fire” data set. For each $q$, the covariates and their
respective penalized coefficient estimates are presented in decreasing order
of the magnitude of the absolute values of the estimators.

As presented in Table 13, it is likely that there are indeed noteworthy
varying covariate effects. Temperature $x_{9}$ appears significant at all
quantile levels but its effects switch sign. The Drought Code $x_{7}$
showcases prominent negative impacts on $y$ from middle to upper-middle $q$,
and the Duff Moisture Code $x_{6}$ appears more influential in positively
affecting $y$ at higher quantile levels. This finding, that higher DMC and
lower DC result in more severe destructions for larger fires, match well to
what [1] found out. We also observed that smaller ($q=0.1$) fires’ impacts
seem to have much fewer key predictors whereas impairments from fires of an
upper-middle scale ($q=0.75$) appear correlated to much more covariates.

###  5.2 Monotone Spline for Quantile Monotonicity Preservation

Consider situations when we have a sole explanatory variable $x$ for a
response variable $y$. Sometimes, we know that $y$ should always decrease or
increase as $x$ increases, yet sample variation may lead to non-monotone
fitted regression curves. At times, fitted curves can present manifest
counterintuitive yet unjustifiable characteristics, casting severe doubt on
the underlying models’ adequacy. One illustrative example is an
aforementioned262626In Section 4.3. generalized gamma model fitted by [32] to
the “WTD” dataset, whose scatter plot is given in Figure 7B. As shown in
Figure 14(a), the GG model’s fitted $0.9^{\text{th}}$ quantile curve notably
shifts upward for large values of $x$, violating the perception that water
table depth should decrease when more rainfall penetrates into the ground
(when flux is bigger). Another case in point is a double kernel nonparametric
smoothing fitting to the aforementionedfootnote 26 “IgG” data set (Figure
21(a)) brought forward by [46]. As in Figure 21(b), [46]’s quantiles exhibit
evident bumps for elder kids, which remain unexplained and disagree with
results obtained by [37] and [24].

(a) Scatter plot

(b) Questionable quantile curves

Figure 21: Scatter plot and fitted quantile curves at
$q=0.05,0.1,0.25,0.5,0.75,0.9,0.95$ by a double kernel method [46] for the
“IgG” data set.

To address such concerns, we opt to impose a monotonicity constraint via
Monotone Cubic Interpolation [14], which ensures monotonicity by the following
two steps:

  1. 1.

Transform $\bm{x}=(x_{1},x_{2},\dots,x_{n})^{T}$ into an $n\times p$ matrix

$\bm{X}_{\operatorname*{bs}}=\begin{pmatrix}{\operatorname*{bs}}_{1}(x_{1})&{\operatorname*{bs}}_{2}(x_{1})&\dots&{\operatorname*{bs}}_{p}(x_{1})\\\
{\operatorname*{bs}}_{1}(x_{2})&{\operatorname*{bs}}_{2}(x_{2})&\dots&{\operatorname*{bs}}_{p}(x_{2})\\\
\vdots&\vdots&\ddots&\vdots\\\
{\operatorname*{bs}}_{1}(x_{n})&{\operatorname*{bs}}_{2}(x_{n})&\dots&{\operatorname*{bs}}_{p}(x_{n})\end{pmatrix}$
(5.11)

using B-spline (bs) with $p-2$ equally spaced knots (2 boundary knots and
$p-4$ inner knots) and intercept = TRUE in R;

  2. 2.

Restrict the regression coefficients
$\bm{\beta}=(\beta_{1},\beta_{2},\dots,\beta_{p})^{T}$ for the model

$\bm{\hat{y}}_{n\times 1}=\bm{X}_{\operatorname*{bs}}\bm{\hat{\beta}}$ (5.12)

to be strictly increasing or decreasing, i.e.,

$\beta_{1}<\beta_{2}<\dots<\beta_{p}\text{ or
}\beta_{1}>\beta_{2}>\dots>\beta_{p}.$ (5.13)

To guarantee the coefficient restriction Equation 5.13 when implementing
monotone regression, we introduce an unrestricted vector of to-be-estimated
regression parameters
$\bm{\theta}=(\theta_{1},\theta_{2},\dots,\theta_{p})^{T}$ specified by

$\displaystyle\theta_{1}=\beta_{1},\theta_{2}=\ln(\beta_{2}-\beta_{1}),\theta_{3}=\ln(\beta_{3}-\beta_{2}),\dots,\theta_{p}=\ln(\beta_{p}-\beta_{p-1})\text{
so that}$ (5.14)
$\displaystyle\beta_{1}=\theta_{1},\beta_{2}=\theta_{1}+e^{\theta_{2}},\beta_{3}=\theta_{1}+e^{\theta_{2}}+e^{\theta_{3}},\dots,\beta_{p}=\theta_{1}+e^{\theta_{2}}+\dots+e^{\theta_{p}}$
(5.15)

for monotonically increasing regression (achieved by requiring
$\beta_{1}<\beta_{2}<\dots<\beta_{p}$) or

$\displaystyle\theta_{1}=\beta_{1},\theta_{2}=\ln(\beta_{1}-\beta_{2}),\theta_{3}=\ln(\beta_{2}-\beta_{3}),\dots,\theta_{p}=\ln(\beta_{p-1}-\beta_{p})\text{
so that}$ (5.16)
$\displaystyle\beta_{1}=\theta_{1},\beta_{2}=\theta_{1}-e^{\theta_{2}},\beta_{3}=\theta_{1}-e^{\theta_{2}}-e^{\theta_{3}},\dots,\beta_{p}=\theta_{1}-e^{\theta_{2}}-\dots-e^{\theta_{p}}$
(5.17)

for monotonically decreasing regression (achieved by requiring
$\beta_{1}>\beta_{2}>\dots>\beta_{p}$).

Getting back to our Quantile Regression MM Algorithm [20] in Section 2.2, we
note that $Q_{\epsilon}(\bm{\theta}|\bm{\theta}^{(t)})$ (Equation 2.7) is not
quadratic in $\bm{\theta}$ anymore since

$f_{i}(\bm{\theta})=\sum_{j=1}^{p}\beta_{j}\cdot{\operatorname*{bs}}_{j}(x_{i})=\begin{cases}\theta_{1}\cdot\sum_{j=1}^{p}{\operatorname*{bs}}_{j}(x_{i})+e^{\theta_{2}}\cdot\sum_{j=2}^{p}{\operatorname*{bs}}_{j}(x_{i})+\dots+e^{\theta_{p}}\cdot{\operatorname*{bs}}_{p}(x_{i})\\\
\theta_{1}\cdot\sum_{j=1}^{p}{\operatorname*{bs}}_{j}(x_{i})-e^{\theta_{2}}\cdot\sum_{j=2}^{p}{\operatorname*{bs}}_{j}(x_{i})-\dots-e^{\theta_{p}}\cdot{\operatorname*{bs}}_{p}(x_{i})\end{cases},$
(5.18)

is not linear in $\bm{\theta}$ anymore. Hence, we can no longer directly solve
for $\bm{\theta}^{(t+1)}$ as in Section 3.1. We resort to a Gauss-Newton
method instead to continuously drive the majorizer downhill. This way, the
descent property (Equation 2.3) still holds and thus the value of our actual
objective function is still guaranteed to decrease iteration after iteration.
Let $r_{i}(\bm{\theta})=y_{i}-f_{i}(\bm{\theta})$ with $f_{i}(\bm{\theta})$ as
specified in Equation 5.18. Then

$\displaystyle\left[\frac{\partial}{\partial\bm{\theta}}r_{i}(\bm{\theta})\right]_{p\times
1}=\left(\frac{\partial
r_{i}(\bm{\theta})}{\partial\theta_{1}},\dots,\frac{\partial
r_{i}(\bm{\theta})}{\partial\theta_{p}}\right)^{T}=\left(\sum_{j=1}^{p}{\operatorname*{bs}}_{j}(x_{i}),\mp
e^{\theta_{2}}\cdot\sum_{j=2}^{p}{\operatorname*{bs}}_{j}(x_{i}),\dots,\mp
e^{\theta_{p}}\cdot{\operatorname*{bs}}_{p}(x_{i})\right)^{T},$
$\displaystyle\left[\frac{\partial^{2}}{\partial\bm{\theta}\partial\bm{\theta}^{T}}r_{i}(\bm{\theta})\right]_{p\times
p}=\begin{pmatrix}\frac{\partial^{2}r_{i}(\bm{\theta})}{\partial\theta_{1}^{2}}&\frac{\partial^{2}r_{i}(\bm{\theta})}{\partial\theta_{2}\theta_{1}}&\dots&\frac{\partial^{2}r_{i}(\bm{\theta})}{\partial\theta_{p}\theta_{1}}\\\
\frac{\partial^{2}r_{i}(\bm{\theta})}{\partial\theta_{1}\theta_{2}}&\frac{\partial^{2}r_{i}(\bm{\theta})}{\partial\theta_{2}^{2}}&\dots&\frac{\partial^{2}r_{i}(\bm{\theta})}{\partial\theta_{p}\theta_{2}}\\\
\vdots&\vdots&\ddots&\vdots\\\
\frac{\partial^{2}r_{i}(\bm{\theta})}{\partial\theta_{1}\theta_{p}}&\frac{\partial^{2}r_{i}(\bm{\theta})}{\partial\theta_{2}\theta_{p}}&\dots&\frac{\partial^{2}r_{i}(\bm{\theta})}{\partial\theta_{p}^{2}}\end{pmatrix}=\begin{pmatrix}0&&&\\\
&e^{\theta_{2}}\cdot\sum_{j=2}^{p}{\operatorname*{bs}}_{j}(x_{i})&&\\\
&&\ddots&\\\
&&&e^{\theta_{p}}\cdot{\operatorname*{bs}}_{p}(x_{i})\end{pmatrix}.$ (5.19)

Denote $\bm{r}_{i}^{\prime}(\bm{\theta})_{p\times 1}$ as
$\frac{\partial}{\partial\bm{\theta}}r_{i}(\bm{\theta})$ and
$\bm{r}_{i}^{\prime\prime}(\bm{\theta})_{p\times p}$ as
$\frac{\partial^{2}}{\partial\bm{\theta}\partial\bm{\theta}^{T}}r_{i}(\bm{\theta})$.
Given current parameter estimate $\bm{\theta}^{(t)}$ and the first, second
differentials

$\displaystyle Q_{\epsilon}^{\prime}(\bm{\theta}|\bm{\theta}^{(t)})_{p\times
1}=$
$\displaystyle\frac{1}{2}\sum_{i=1}^{n}\left[\frac{r_{i}(\bm{\theta})}{\epsilon+|r_{i}(\bm{\theta}^{(t)})|}+2q-1\right]\bm{r}_{i}^{\prime}(\bm{\theta})\text{
and}$ $\displaystyle
Q_{\epsilon}^{\prime\prime}(\bm{\theta}|\bm{\theta}^{(t)})_{p\times p}=$
$\displaystyle\frac{1}{2}\sum_{i=1}^{n}\left[\frac{1}{\epsilon+|r_{i}(\bm{\theta}^{(t)})|}\bm{r}_{i}^{\prime}(\bm{\theta})\bm{r}_{i}^{\prime}(\bm{\theta})^{T}+\left(\frac{r_{i}(\bm{\theta})}{\epsilon+|r_{i}(\bm{\theta}^{(t)})|}+2q-1\right)\bm{r}_{i}^{\prime\prime}(\bm{\theta})\right],$
(5.20)

we can calculate the Gauss-Newton step direction at iteration t as the
$p\times 1$ vector

$\bm{\Delta}_{\epsilon}^{(t)}=-\left[\bm{r}^{\prime}\big{(}\bm{\theta}^{(t)}\big{)}^{T}\bm{W}_{\epsilon}\big{(}\bm{\theta}^{(t)}\big{)}\bm{r}^{\prime}\big{(}\bm{\theta}^{(t)}\big{)}+\bm{H}\big{(}\bm{\theta}^{(t)}\big{)}\right]^{-1}\bm{r}^{\prime}\big{(}\bm{\theta}^{(t)}\big{)}^{T}\bm{v}_{\epsilon}\big{(}\bm{\theta}^{(t)}\big{)}^{T},$
(5.21)

where

$\displaystyle\bm{r}^{\prime}(\bm{\theta})_{n\times
p}=\begin{pmatrix}\bm{r}^{\prime}_{1}(\bm{\theta})^{T}\\\
\bm{r}^{\prime}_{2}(\bm{\theta})^{T}\\\ \vdots\\\
\bm{r}^{\prime}_{n}(\bm{\theta})^{T}\\\
\end{pmatrix},\quad\bm{W}_{\epsilon}\big{(}\bm{\theta}^{(t)}\big{)}_{n\times
n}=\begin{pmatrix}\frac{1}{\epsilon+|r_{1}(\bm{\theta}^{(t)})|}&&&\\\
&\frac{1}{\epsilon+|r_{2}(\bm{\theta}^{(t)})|}&&\\\ &&\ddots&\\\
&&&\frac{1}{\epsilon+|r_{n}(\bm{\theta}^{(t)})|}\end{pmatrix},$
$\displaystyle\bm{H}\big{(}\bm{\theta}|\bm{\theta}^{(t)}\big{)}_{p\times
p}=\sum_{i=1}^{n}\left(\frac{r_{i}(\bm{\theta})}{\epsilon+|r_{i}(\bm{\theta}^{(t)})|}+2q-1\right)\bm{r}_{i}^{\prime\prime}(\bm{\theta}),\text{
and}$
$\displaystyle\bm{v}_{\epsilon}\big{(}\bm{\theta}|\bm{\theta}^{(t)}\big{)}_{1\times
n}=\left(\frac{r_{1}(\bm{\theta})}{\epsilon+|r_{1}(\bm{\theta}^{(t)})|}+2q-1,\dots,\frac{r_{n}(\bm{\theta})}{\epsilon+|r_{n}(\bm{\theta}^{(t)})|}+2q-1\right).$
(5.22)

The next parameter estimate given by

$\displaystyle\bm{\theta}^{(t+1)}=\bm{\theta}^{(t)}+\alpha^{(t)}\bm{\Delta}_{\epsilon}^{(t)},\text{
where}$ $\displaystyle\alpha^{(t)}=\begin{cases}1,\;\text{ if
}\bm{\Delta}_{\epsilon}^{(t)}=\bm{0}_{p\times 1},\text{ i.e., a stationary
point}\\\
\max\\{2^{-v}:Q_{\epsilon}(\bm{\theta}^{(t)}+2^{-v}\bm{\Delta}_{\epsilon}^{(t)}|\bm{\theta}^{(t)})<Q_{\epsilon}(\bm{\theta}^{(t)}|\bm{\theta}^{(t)}),v\in\mathbb{N}\\},\text{
otherwise}\end{cases}$ (5.23)

then ensures $Q_{\epsilon}(\bm{\theta}^{(t+1)}|\bm{\theta}^{(t)})\leq
Q_{\epsilon}(\bm{\theta}^{(t)}|\bm{\theta}^{(t)})$.272727When implementing the
Gauss-Newton approach in statistical software like R, we take $\alpha^{(t)}=1$
(interpreted as exploring some new terrain) if
$Q_{\epsilon}(\bm{\theta}^{(t)}+2^{-v}\bm{\Delta}_{\epsilon}^{(t)}|\bm{\theta}^{(t)})\geq
Q_{\epsilon}(\bm{\theta}^{(t)}|\bm{\theta}^{(t)})\;\forall\;v\in\mathbb{N},v\leq
M$ for some suitably big $M\in\mathbb{N}$.

We implemented this MM monotone quantile regression algorithm on the “WTD” and
“IgG” data sets mentioned at the beginning of this subsection. As Figure 22
suggests, monotone spline interpolation achieved via MM successfully preserved
our desired monotonicity while capturing the data sets’ patterns reasonably
well.

Figure 22: Quantile curves (green lines) fitted using monotone B-spline with
one inner knot at $q=0.9$ for the “WTD” dataset (A) and at
$q=0.1,0.25,0.5,0.75,0.9$ for the “IgG” dataset (B), respectively. The purple
lines are the corresponding rq empirical estimates.

##  6 Discussion

This article explores MM algorithms in diverse quantile regression settings.
Applications of MM to several real data sets and two sets of simulation
studies comparing MM with existing tested methods showcase our MM algorithms’
effectiveness.

One main advance we have made is putting forward a novel parametric fitting
method that efficiently models interested quantiles simultaneously for
nonparametric quantile regression, which is achieved via MM detailed in
Section 4.1. MM’s computation speed edge over some other prominent algorithms
when estimating coefficients at different quantiles separately [35] has also
been supported by our simulation study in Section 3.3.

One major flaw of our MM is that it gets relatively slow compared to qrcm when
estimating different quantile coefficients simultaneously, which should result
from having to invert an $hp\times hp$ matrix
$\sum_{a=1}^{k}\bm{D}^{T}(q_{a})\bm{W}_{a}^{(t)}\bm{D}(q_{a})$ (in Section
3.2) at each iteration. It is, however, quite promising to address this
current weakness, as some authors have successfully invented fast MM
algorithms that bypass large matrix inversion problems ([29]; [12]). In
situations where super fast algorithms are needed, a finite smoothing
algorithm raised by [6], which operates by directly minimizing a nice smooth
approximate function of
$L_{q}(\bm{\theta})=\sum_{i=1}^{n}\rho_{q}(r_{i}(\bm{\theta}))$, may also come
to the rescue.

##  Appendix A Additional Proofs

Proposition A.1.

$\displaystyle
Q(\bm{\theta}|\bm{\theta}^{(t)})=\sum_{i=1}^{n}\zeta_{q}(r_{i}(\bm{\theta})|r_{i}(\bm{\theta}^{(t)}))=\sum_{i=1}^{n}\frac{1}{4}\left\\{\frac{r_{i}^{2}(\bm{\theta})}{|r_{i}(\bm{\theta}^{(t)})|}+(4q-2)r_{i}(\bm{\theta})+|r_{i}(\bm{\theta}^{(t)})|\right\\}\text{
majorizes }$ $\displaystyle
L(\bm{\theta})=\sum_{i=1}^{n}\rho_{q}(r_{i}(\bm{\theta}))=\sum_{i=1}^{n}qr_{i}(\bm{\theta})-r_{i}(\bm{\theta})\mathbbm{1}_{\\{r_{i}(\bm{\theta})<0\\}}\text{
at }\bm{\theta}=\bm{\theta}^{(t)}.$

###### Proof.

$\displaystyle
Q(\bm{\theta}^{(t)}|\bm{\theta}^{(t)})-L(\bm{\theta}^{(t)})=\sum_{i=1}^{n}\left[\zeta_{q}(r_{i}(\bm{\theta}^{(t)})|r_{i}(\bm{\theta}^{(t)}))-\rho_{q}(r_{i}(\bm{\theta}^{(t)}))\right]$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{n}\frac{1}{4}\left[|r_{i}(\bm{\theta}^{(t)})|-2r_{i}(\bm{\theta}^{(t)})+|r_{i}(\bm{\theta}^{(t)})|\right]+r_{i}(\bm{\theta}^{(t)})\mathbbm{1}_{\\{r_{i}(\bm{\theta}^{(t)})<0\\}}$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{n}\left\\{r_{i}(\bm{\theta}^{(t)})\left(-\frac{1}{2}+\mathbbm{1}_{\\{r_{i}(\bm{\theta}^{(t)})<0\\}}\right)+\frac{1}{2}|r_{i}(\bm{\theta}^{(t)})|\right\\}=0.$

It remains to show $Q(\bm{\theta}|\bm{\theta}^{(t)})$ satisfies Equation 2.1.
$\forall\text{ residual }r,$ letting

$\displaystyle h(r)=\,$
$\displaystyle\zeta_{q}(r|r^{(t)})-\rho_{q}(r)=\frac{1}{4}\left(\frac{r^{2}}{|r^{(t)}|}-2r+|r^{(t)}|\right)+r\mathbbm{1}_{\\{r<0\\}},\text{
we have }$ $\displaystyle h^{\prime}(r)=\,$
$\displaystyle\zeta_{q}^{\prime}(r|r^{(t)})-{\rho_{q}}^{\prime}(r)=\frac{1}{2}\left(\frac{r}{|r^{(t)}|}-1\right)+\mathbbm{1}_{\\{r<0\\}}=\frac{1}{2}\left(\frac{r-|r^{(t)}|(1-2\cdot\mathbbm{1}_{\\{r<0\\}})}{|r^{(t)}|}\right)$
$\displaystyle=\,$
$\displaystyle\begin{cases}\frac{1}{2}\cdot\frac{r-|r^{(t)}|}{|r^{(t)}|}\;\text{
if }r\geq 0\\\ \frac{1}{2}\cdot\frac{r+|r^{(t)}|}{|r^{(t)}|}\;\text{ if
}r<0,\qquad\text{and}\\\ \end{cases}$ $\displaystyle h^{\prime\prime}(r)=$
$\displaystyle\frac{1}{2|r^{(t)}|}>0\;(r^{(t)}\neq 0)\quad\forall\;r.$ (A.1)

Observe from Appendix A that $h(\pm r^{(t)})=h^{\prime}(\pm r^{(t)})=0$ and

$\displaystyle h^{\prime}(r)\leq 0\>\text{ if
}r\in(-\infty,-|r^{(t)}|]\cup[0,|r^{(t)}|]\text{ and }\geq 0\>\text{ if
}r\in[-|r^{(t)}|,0]\cup[|r^{(t)}|,\infty)$ $\displaystyle\iff$ $\displaystyle
h(r)\>\text{ is non-increasing for }r\in(-\infty,-|r^{(t)}|]\cup[0,|r^{(t)}|]$
$\displaystyle\text{ and is non-decreasing for
}r\in[-|r^{(t)}|,0]\cup[|r^{(t)}|,\infty)$ $\displaystyle\Rightarrow\>$
$\displaystyle h(r)\geq\min\\{h(r^{(t)}),h(-r^{(t)})\\}=0\;\forall\>r$
$\displaystyle\iff$
$\displaystyle\zeta_{q}(r|r^{(t)})\geq\rho_{q}(r)\;\forall\>r.$

Hence $Q(\bm{\theta|\theta}^{(t)})\geq L(\bm{\theta})\;\forall\;\bm{\theta}$
and Equation 2.1 is satisfied. ∎

Proposition A.2. Consider the linear quantile regression model

$Y_{i}=\bm{X}_{i}^{T}\bm{\beta}_{q}+U_{i},\>i=1,2,\dots,n,\>q\in(0,1)$ (A.2)

with loss $\rho_{q}(u)=|u|[q\mathbbm{1}_{\\{u\geq
0\\}}+(1-q)\mathbbm{1}_{\\{u<0\\}}]=qu-u\mathbbm{1}_{\\{u<0\\}}$ and error
distribution

$U_{1},\dots,U_{n}\overset{\operatorname*{iid}}{\sim}\text{ an asymmetric
Laplace distribution with PDF
}f_{U}(u)=\frac{q(1-q)}{\sigma}e^{-\frac{\rho_{q}(u)}{\sigma}}.$ (A.3)

Given a fixed $\lambda\geq 0$ and an obtained (penalized) coefficient
estimator $\hat{\bm{\beta}}_{q}=\hat{\bm{\beta}}_{\lambda}(q)$ that attempts
to minimize $L_{q}^{\epsilon}(\bm{\beta})$ in Equation 5.4, a BIC is given by

${\operatorname*{BIC}}_{q}\left(\hat{\bm{\beta}}_{q}\right)=\ln\left(\sum_{i=1}^{n}\rho_{q}\left(r_{i}\left(\hat{\bm{\beta}}_{q}\right)\right)\right)+|\mathit{S}_{q}|\frac{\ln(n)}{2n}$
(A.4)

if we approximate $\sigma$ by $\hat{\sigma}_{\operatorname*{mle}}$.

###### Proof.

Given sample observations $\\{(\bm{x}_{i},y_{i}):1\leq i\leq n\\}$, the log-
likelihood is

$l(\bm{\beta}_{q},\sigma)=\ln\left[\prod_{i=1}^{n}f_{U}(y_{i}-\bm{x}_{i}^{T}\bm{\beta}_{q})\right]=-n\ln(\sigma)-\sum_{i=1}^{n}\frac{\rho_{q}(y_{i}-\bm{x}_{i}^{T}\bm{\beta}_{q})}{\sigma}+\operatorname*{const}.$
(A.5)

For any given $\hat{\bm{\beta}}_{q}$, we let
$g(\sigma)=l\big{(}\hat{\bm{\beta}}_{q},\sigma\big{)},\>\sigma\in(0,\infty)$.
Since $g$ is differentiable everywhere on $(0,\infty)$ and

$\displaystyle g^{\prime}(\sigma)$ $\displaystyle=\frac{\partial
l\left(\bm{\beta}_{q},\sigma\right)}{\partial\sigma}\bigg{|}_{\bm{\beta}_{q}=\hat{\bm{\beta}}_{q}}=-\frac{n}{\sigma}+\sum_{i=1}^{n}\frac{\rho_{q}\left(y_{i}-\bm{x}_{i}^{T}\hat{\bm{\beta}}_{q}\right)}{\sigma^{2}}=\frac{-n\sigma+\sum_{i=1}^{n}\rho_{q}\left(y_{i}-\bm{x}_{i}^{T}\hat{\bm{\beta}}_{q}\right)}{\sigma^{2}}$
has a unique zero (and hence the unique critical point by Fermat’s Theorem) at
$\displaystyle\sigma$
$\displaystyle=\hat{\sigma}\big{(}\hat{\bm{\beta}_{q}}\big{)}=\frac{\sum_{i=1}^{n}\rho_{q}\left(y_{i}-\bm{x}_{i}^{T}\hat{\bm{\beta}}_{q}\right)}{n}$
(A.6)

with $g^{\prime}(\sigma)>0$ for
$\sigma<\hat{\sigma}\big{(}\hat{\bm{\beta}}_{q}\big{)}$ and
$g^{\prime}(\sigma)<0$ for
$\sigma>\hat{\sigma}\big{(}\hat{\bm{\beta}_{q}}\big{)}$,
$\sigma=\hat{\sigma}\big{(}\hat{\bm{\beta}}_{q}\big{)}$ given in Appendix A is
the unique local maxima. Since $g(\sigma)\rightarrow-\infty$ as
$\sigma\rightarrow\infty$ or as $\sigma\rightarrow 0_{+}$,

$\hat{\sigma}_{\operatorname*{mle}}\left(\hat{\bm{\beta}}_{q}\right)=\frac{\sum_{i=1}^{n}\rho_{q}\left(y_{i}-\bm{x}_{i}^{T}\hat{\bm{\beta}}_{q}\right)}{n}\in(0,\infty)\text{
satisfies
}g\left(\hat{\sigma}_{\operatorname*{mle}}\left(\hat{\bm{\beta}}_{q}\right)\right)=\underset{\sigma\in(0,\infty)}{\max}\,g(\sigma)=\underset{\sigma\in(0,\infty)}{\max}\,l\left(\hat{\bm{\beta}}_{q},\sigma\right)$

and is the unique MLE of $\sigma$ given $\hat{\bm{\beta}}_{q}$. Setting
$\sigma=\hat{\sigma}_{\operatorname*{mle}}\left(\hat{\bm{\beta}}_{q}\right)=\frac{\sum_{i=1}^{n}\rho_{q}\left(y_{i}-\bm{x}_{i}^{T}\hat{\bm{\beta}}_{q}\right)}{n}$,
we get282828note that $n$ is given and hence constant

$\displaystyle{\operatorname*{BIC}}_{q}(\hat{\bm{\beta}}_{q})$
$\displaystyle=-2l\left(\hat{\bm{\beta}}_{q},\hat{\sigma}_{\operatorname*{mle}}(\hat{\bm{\beta}}_{q})\right)+|\mathit{S}_{q}|\ln(n)=2n\ln\left[\hat{\sigma}_{\operatorname*{mle}}(\hat{\bm{\beta}}_{q})\right]+2n-2\operatorname*{const}+|\mathit{S}_{q}|\ln(n)$
$\displaystyle=2n\left(\ln\left[\frac{\sum_{i=1}^{n}\rho_{q}(y_{i}-\bm{x}_{i}^{T}\hat{\bm{\beta}}_{q})}{n}\right]+|\mathit{S}_{q}|\frac{\ln(n)}{2n}\right)+2n-2\operatorname*{const}$
$\displaystyle=2n\left(\ln\left[\sum_{i=1}^{n}\rho_{q}(y_{i}-\bm{x}_{i}^{T}\hat{\bm{\beta}}_{q})\right]+|\mathit{S}_{q}|\frac{\ln(n)}{2n}\right)+2n-2\operatorname*{const}-2n\ln(n)$
$\displaystyle={\operatorname*{const}}_{1}\left(\ln\left[\sum_{i=1}^{n}\rho_{q}(y_{i}-\bm{x}_{i}^{T}\hat{\bm{\beta}}_{q})\right]+|\mathit{S}_{q}|\frac{\ln(n)}{2n}\right)+{\operatorname*{const}}_{2},\text{
where }{\operatorname*{const}}_{1}=2n>0.$ (A.7)

Hence, Appendix A is equivalent to Equation A.4 as a BIC metric. ∎

##  Appendix B Additional Tables and Figures

Figure B.1: Estimated regression coefficients $\hat{\bm{\beta}}(q)$ versus quantile levels $q$ plot for the “pollution” data set and our model with the logistic basis. The orange lines are rq’s empirical regression coefficient estimates at 1000 quantile levels. The blue and green lines denote estimated $\hat{\bm{\beta}}(q)$ produced by qrcm and our MM algorithm described in Section 3.2, respectively. Figure B.2: Estimated regression coefficients $\hat{\bm{\beta}}(q)$ versus quantile levels $q$ plot for the “pollution” data set and our model with natural spline basis (with knots $0.25,0.5,0.75$). The orange lines are rq’s empirical regression coefficient estimates at 1000 quantile levels. The blue and green lines denote estimated $\hat{\bm{\beta}}(q)$ produced by qrcm and our MM algorithm described in Section 3.2, respectively. Figure B.3: Estimated regression coefficients $\hat{\bm{\beta}}(q)$ versus quantile levels $q$ plot for the “pollution” data set and our model with natural spline basis (with knots $0.1,0.5,0.9$). The orange lines are rq’s empirical regression coefficient estimates at 1000 quantile levels. The blue and green lines denote estimated $\hat{\bm{\beta}}(q)$ produced by qrcm and our MM algorithm described in Section 3.2, respectively. Figure B.4: Estimated regression coefficients $\hat{\bm{\beta}}(q)$ versus quantile levels $q$ plot for the “pollution” data set and our model with natural spline basis (with knots $0.1,0.4,0.6,0.9$). The orange lines are rq’s empirical regression coefficient estimates at 1000 quantile levels. The blue and green lines denote estimated $\hat{\bm{\beta}}(q)$ produced by qrcm and our MM algorithm described in Section 3.2, respectively. Figure B.5: Estimated regression coefficients $\hat{\bm{\beta}}(q)$ versus quantile levels $q$ plot for the “pollution” data set and our model with natural spline basis (with knots $0.1,0.3,0.7,0.9$). The orange lines are rq’s empirical regression coefficient estimates at 1000 quantile levels. The blue and green lines denote estimated $\hat{\bm{\beta}}(q)$ produced by qrcm and our MM algorithm described in Section 3.2, respectively. Figure B.6: Estimated regression coefficients $\hat{\bm{\beta}}(q)$ versus quantile levels $q$ plot for the “pollution” data set and our model with natural spline basis (with knots $0.2,0.4,0.6,0.8$). The orange lines are rq’s empirical regression coefficient estimates at 1000 quantile levels. The blue and green lines denote estimated $\hat{\bm{\beta}}(q)$ produced by qrcm and our MM algorithm described in Section 3.2, respectively. Figure B.7: Estimated regression coefficients $\hat{\bm{\beta}}(q)$ versus quantile levels $q$ plot for the “pollution” data set and our model with natural spline basis (with knots $0.1,0.3,0.5,0.7,0.9$). The orange lines are rq’s empirical regression coefficient estimates at 1000 quantile levels. The blue and green lines denote estimated $\hat{\bm{\beta}}(q)$ produced by qrcm and our MM algorithm described in Section 3.2, respectively. |  |  |  |  |
---|---|---|---|---|---
|  |  |  |  |
|  |  |  |  |
|  |  |  |  |
|  |  |  |  |
|  |  |  |  |
|  |  |  |  |
|  |  |  |  |
|  |  |  |  |
|  |  |  |  |
Figure B.8: Conditional $\hat{y}$ quantile plots given the 60 sample observations in the “pollution” data set [31] – logistic basis. The orange lines are rq’s empirical predictions at $1000$ quantile levels; the green lines represent MM’s conditional quantile prediction functions $\hat{Q}_{q}(Y|\bm{x})=\hat{y}(q|\bm{x})$, $q\in(0,1)$; the blue lines denote qrcm’s predicted values $\hat{y}(q|\bm{x})$. |  |  |  |  |
---|---|---|---|---|---
|  |  |  |  |
|  |  |  |  |
|  |  |  |  |
|  |  |  |  |
|  |  |  |  |
|  |  |  |  |
|  |  |  |  |
|  |  |  |  |
|  |  |  |  |
Figure B.9: Conditional $\hat{y}$ quantile plots given the 60 sample
observations in the “pollution” data set [31] – natural spline basis with
knots $0.2,0.4,0.6,0.8$. The orange lines are rq’s empirical predictions at
$1000$ quantile levels; the green lines represent MM’s conditional quantile
prediction functions $\hat{Q}_{q}(Y|\bm{x})=\hat{y}(q|\bm{x})$, $q\in(0,1)$;
the blue lines denote qrcm’s predicted values $\hat{y}(q|\bm{x})$. Figure
B.10: Fitted quantile curves for the “triceps” data set - natural spline
transformation on $x$ with 3 equally spaced knots. We have adopted the
logistic basis (A) and natural spline bases with 3 (B), 4 (C), 5 (D), and 6
(E) equally spaced knots for quantile levels. Figure B.11: Fitted quantile
curves for the “triceps” data set - natural spline transformation on $x$ with
4 equally spaced knots. We have adopted the logistic basis (A) and natural
spline bases with 3 (B), 4 (C), 5 (D), and 6 (E) equally spaced knots for
quantile levels. Figure B.12: Fitted quantile curves for the “triceps” data
set - natural spline transformation on $x$ with 5 equally spaced knots. We
have adopted the logistic basis (A) and natural spline bases with 3 (B), 4
(C), 5 (D), and 6 (E) equally spaced knots for quantile levels. Figure B.13:
Fitted quantile curves for the “triceps” data set - natural spline
transformation on $x$ with 7 equally spaced knots. We have adopted the
logistic basis (A) and natural spline bases with 3 (B), 4 (C), 5 (D), and 6
(E) equally spaced knots for quantile levels. Figure B.14: Fitted quantile
curves for the “triceps” data set - natural spline transformation on $x$ with
8 equally spaced knots. We have adopted the logistic basis (A) and natural
spline bases with 3 (B), 4 (C), 5 (D), and 6 (E) equally spaced knots for
quantile levels. Figure B.15: Fitted quantile curves for the “WTD” data set -
natural spline transformation on $x$ with 3 equally spaced knots. We have
adopted the logistic basis (A) and natural spline bases with 3 (B), 4 (C), 5
(D), and 6 (E) equally spaced knots for quantile levels. Figure B.16: Fitted
quantile curves for the “WTD” data set - natural spline transformation on $x$
with 4 equally spaced knots. We have adopted the logistic basis (A) and
natural spline bases with 3 (B), 4 (C), 5 (D), and 6 (E) equally spaced knots
for quantile levels. Figure B.17: Fitted quantile curves for the “WTD” data
set - natural spline transformation on $x$ with 5 equally spaced knots. We
have adopted the logistic basis (A) and natural spline bases with 3 (B), 4
(C), 5 (D), and 6 (E) equally spaced knots for quantile levels. Figure B.18:
Fitted quantile curves for the “WTD” data set - natural spline transformation
on $x$ with 6 equally spaced knots. We have adopted the logistic basis (A) and
natural spline bases with 3 (B), 4 (C), 5 (D), and 6 (E) equally spaced knots
for quantile levels. Figure B.19: Fitted quantile curves for the “WTD” data
set - natural spline transformation on $x$ with 8 equally spaced knots. We
have adopted the logistic basis (A) and natural spline bases with 3 (B), 4
(C), 5 (D), and 6 (E) equally spaced knots for quantile levels.

Quantile Level | 0.1 | 0.25 | 0.5 | 0.75 | 0.9
---|---|---|---|---|---
MM xseq3 | 0.00705 | 0.00548 | 0.00596 | 0.00768 | 0.01510
MM xseq4 | 0.00850 | 0.00671 | 0.00759 | 0.00936 | 0.01754
MM xseq5 | 0.00970 | 0.00782 | 0.00895 | 0.01066 | 0.01973
MM xseq6 | 0.01133 | 0.00877 | 0.01002 | 0.01198 | 0.02253
MM xseq7 | 0.01278 | 0.00974 | 0.01116 | 0.01341 | 0.02527
MM xasym7 | 0.00783 | 0.00569 | 0.00611 | 0.00793 | 0.01536
MM xseq3qseq3 | 0.00825 | 0.00589 | 0.00666 | 0.00984 | 0.01744
MM xseq3qseq5 | 0.00680 | 0.00621 | 0.00696 | 0.00818 | 0.01566