#  SPQR: An R Package for Semi-Parametric Density and Quantile Regression

by Steven G. Xu, Reetam Majumder and Brian J. Reich

###### Abstract

We develop an R package SPQR that implements the semi-parametric quantile
regression (SPQR) method in Xu and Reich (2021). The method begins by fitting
a flexible density regression model using monotonic splines whose weights are
modeled as data-dependent functions using artificial neural networks.
Subsequently, estimates of conditional density and quantile process can all be
obtained. Unlike many approaches to quantile regression that assume a linear
model, SPQR allows for virtually any relationship between the covariates and
the response distribution including non-linear effects and different effects
on different quantile levels. To increase the interpretability and
transparency of SPQR, model-agnostic statistics developed by Apley and Zhu
(2020) are used to estimate and visualize the covariate effects and their
relative importance on the quantile function. In this article, we detail how
this framework is implemented in SPQR and illustrate how this package should
be used in practice through simulated and real data examples.

## Introduction

Quantile regression (QR) (Koenker and Bassett Jr, 1978) is a widely used tool
in applications where covariate effect on non-central part of the response
distribution is of interest. Such problems routinely arrive in economics,
health and environment studies where the response is often heavy-tailed and/or
heteroscatistic. Although modeling a single quantile or a small subset of
quantiles often serves as the primary interest in most applications, the full
potential of QR lies in estimation of the full quantile process to obtain a
comprehensive description of the varying covariate effects across the quantile
domain.

Most existing implementations of quantile regression in R assume a linear
relationship between the conditional quantile and covariates, examples include
bayesQR (Benoit and Van den Poel, 2017), lqr (Galarza et al., 2022) and
quantreg Koenker (2022). These methods become too restrictive, leading to
large estimation bias, when the underlying relationship is highly nonlinear.
To overcome this challenge, semi-parametric QR models that allow some or some
transformation of the covariate effects be modeled flexibly have been
proposed, examples are quantreg.nonpar (Lipsitz et al., 2016), plaqr (Maidman,
2017) and qgam (Fasiolo et al., 2021) that use generalized additive models;
and siqr (Yu, 2021) that uses a single index model. Additive models allow
flexible and interpretable modeling of main effects, but explicitly specifying
the interaction terms quickly become a tedious task as the data dimension
increases. Single index model alleviates curse of dimensionality by projecting
the covariates to a 1-dimensional space, but when the intrinsic data dimension
is high such projection can be too coarse to retain all valuable information.
A more serious disadvantage shared by all aforementioned methods is that they
do not ensure the monotonicity of the estimated conditional quantile function.
Without such constraint, the estimated quantiles can cross each other when the
sample size is small, leading to a invalid response distribution that is
difficult or impossible to draw inference from.

Few existing packages are capable of estimating multiple quantiles under non-
crossing constraints. The qrjoint package (Tokdar and Cunningham, 2019)
directly models the linear coefficient function as a monotone process and
allows simultaneous estimation of non-crossing quantile planes. The mcqrnn()
function in the qrnn package (Cannon, 2018) can model non-linear non-crossing
quantile curves by treating the quantile level as an observed covariate and
imposing positivity constraints on its correspond weights in the neural
network. However, the model suffers from high variance when the sample size is
small (Xu and Reich, 2021) and uncertainty quantification is not
straightforward.

The purpose of this article is to introduce a new R package, named SPQR (Xu
and Majumder, 2022), for flexible QR modeling. This package implements the
semi-parametric quantile regression (SPQR) model proposed by Xu and Reich
(2021) which allows simultaneous estimation of non-linear non-crossing
quantile curves. The method begins by specifying a semi-parametric model for
the conditional distribution function using shape-constrained splines, the
coefficients are then modeled non-parametrically as functions of covariates
using artificial neural networks. As a result, valid estimates of the
conditional response distribution and its quantiles can be simultaneously
obtained. The SPQR provides three approaches for fitting SPQR: maximum
likelihood estimation (MLE) and maximum _a posteriori_ probability (MAP) which
provide point estimates, and Markov chain Monte Carlo (MCMC) which provides
uncertainty quantification through posterior samples. To our best knowledge
and at the time of writing, SPQR is the only package that implements a fully
Bayesian framework for estimating semi-parametric non-crossing QR. The main
computations of SPQR rely on torch (Falbel and Luraschi, 2022) and Rcpp
(Eddelbuettel et al., 2022a). Specifically, MLE and MAP are optimized using
the Adam routine (Kingma and Ba, 2014) in torch and take advantage of its GPU
capability, and MCMC is implemented using Rcpp and RcppArmadillo (Eddelbuettel
et al., 2022b) for efficient gradient computation. To increase the
interpretability and transparency of SPQR, SPQR also provides function to
compute quantile accumulative local effects (ALE; Apley and Zhu, 2020) to
characterize and visualize quantile-dependent main and interaction effects, as
well as function to compute quantile-dependent variable importance measures.

The rest of the article is organized as follows. We first give a brief review
of the methodology background of SPQR and quantile ALEs. We then introduce the
implemented computation approaches for fitting SPQR. The subsequent section
introduces the SPQR package and its main features in detail. Finally, we
demonstrate the usage and effectiveness of SPQR in density and quantile
regression problems using simulated and real data. In particular, we apply our
package to the Australia electric demand data, available in the qgam package,
and analyze the quantile effects of temperature and time on residential
electric energy consumption. The last section concludes the paper.

## Methodology background

### Density regression model

To estimate quantile effects, we first build a model for the probability
density function (PDF) of the response $Y$ conditioned on covariates
$\mbox{\bf X}=(X_{1},...,X_{p})$. We assume that the response and all $p$
covariates are scaled to $[0,1]$ to simplify the specifications of prior
distributions and basis functions. The PDF is assumed to be a function of $K$
second-order M-spline basis functions, $M_{1}(y),...,M_{K}(y)$, with equally-
spaced knots spanning $[0,1]$. Each M-spline basis function is itself a valid
PDF on $[0,1]$ (Ramsay, 1988), and therefore any convex combination of these
functions is also a valid PDF. The model is

$f(y|\mbox{\bf X})=\sum_{k=1}^{K}\theta_{k}(\mbox{\bf X})M_{k}(y),$ (1)

where the probabilities $\theta_{k}(\mbox{\bf X})$ satisfy
$\theta_{k}(\mbox{\bf X})\geq 0$ and $\sum_{k=1}^{K}\theta_{k}(\mbox{\bf
X})=1$ for all possible X.

To provide flexibility and ensure non-negative weights $\theta_{k}(\mbox{\bf
X})$ that sum to one, we use a fully connected neural network (NN) with
softmax output activation. An NN with $L$ layers ($L-1$ hidden layers) is
parameterized by a set of weight matrices $\mathcal{W}=\\{\mbox{\bf
W}^{(1)},...,\mbox{\bf W}^{(L)}\\}$, with $\mbox{\bf
W}^{(l)}\in\mathbb{R}^{(V_{l}+1)\times V_{l-1}}$ where $V_{l}$ is the number
of units (excluding the intercept/bias node) in layer $l$. We define the input
layer as layer $0$ for simplicity. The first layer is

$z^{(1)}_{i}(\mbox{\bf
X},\mathcal{W})=W^{(1)}_{i0}+\sum_{j=1}^{p}W^{(1)}_{ij}X_{j}.$ (2)

We make the functional dependence on X and $\mathcal{W}$ explicit in our
notation as it will help clarify what follows. Subsequent hidden layers are
defined by the recursion

$\begin{split}u_{i}^{(l)}(\mbox{\bf
X},\mathcal{W})&=\phi\left\\{z^{(l)}_{i}(\mbox{\bf
X},\mathcal{W})\right\\}\mbox{\ \ \ for\ \ }l\in\\{1,...,L-1\\}\\\
z^{(l+1)}_{i}(\mbox{\bf
X},\mathcal{W})&=W_{i0}^{(l+1)}+\sum_{j=1}^{V_{l}}W^{(l+1)}_{ij}u_{j}^{(l)}(\mbox{\bf
X},\mathcal{W})\end{split}$ (3)

where $\phi$ is the activation function taken to be either the hyperbolic
tangent function $\phi(u)=(e^{2u}-1)/(e^{2u}+1)$ or the rectified-linear
function $\phi(u)=\max(0,u)$. Finally, the softmax activation is used in the
output layer to ensure probabilities sum to one

$\theta_{k}(\mbox{\bf X},\mathcal{W})=\frac{\exp\left\\{z^{(L)}_{k}(\mbox{\bf
X},\mathcal{W})\right\\}}{\sum_{i=1}^{K}\exp\left\\{z^{(L)}_{i}(\mbox{\bf
X},\mathcal{W})\right\\}}.$ (4)

The SPQR model for the conditional PDF is then

$f(y|\mathcal{W},\mbox{\bf
X})=\sum_{k=1}^{K}M_{k}(y)\frac{\exp\left\\{z^{(L)}_{k}(\mbox{\bf
X},\mathcal{W})\right\\}}{\sum_{i=1}^{K}\exp\left\\{z^{(L)}_{i}(\mbox{\bf
X},\mathcal{W})\right\\}}.$ (5)

Endowed with the approximation theories of spline and NN, the model in (5) can
approximate any smooth PDF as $K$ and the $V_{l}$ increase. Therefore, this
model can capture complex relationships, such as covariate-dependent variance,
skewness, or likelihood of extreme events.

Computing the PDF is simple and fast given the parameters $\mathcal{W}$, and
since the integral of M-spline functions are I-spline functions (Ramsay,
1988), (5) immediately gives rise to an expression for the cumulative
distribution function (CDF)

$F(y|\mathcal{W},\mbox{\bf
X})=\sum_{k=1}^{K}I_{k}(y)\frac{\exp\left\\{z^{(L)}_{k}(\mbox{\bf
X},\mathcal{W})\right\\}}{\sum_{i=1}^{K}\exp\left\\{z^{(L)}_{i}(\mbox{\bf
X},\mathcal{W})\right\\}},$ (6)

where $I_{k}(y)$ are I-spline basis functions. The conditional quantile
function for quantile level $\tau\in(0,1)$ is defined as the function
$Q(\tau|\mathcal{W},\mbox{\bf X})$ so that

$F\\{Q(\tau|\mathcal{W},\mbox{\bf X})|\mathcal{W},\mbox{\bf X}\\}=\tau.$

The conditional quantile function for this model is not available in closed-
form, but can be approximated by numerically inverting
$F(y|\mathcal{W},\mbox{\bf X})$. Given that (6) models a valid CDF, the
conditional quantile function estimated through this approach satisfies the
non-crossing constraint

$\frac{\partial Q(\tau|\mathcal{W},\mbox{\bf X})}{\partial\tau}>0,\ \forall\
\mbox{\bf X}$

and does not require any second-stage monotonization treatment.

### Summarizing covariate effects on quantiles

SPQR has the advantage of being a flexible semi-parametric model that can
capture complex non-linear covariate effects on various aspects of the
response distribution. A disadvantage is that it is difficult to interpret
individual parameters because the weights $\mathcal{W}$ are not individually
identified and do not correspond to meaningful quantities. In most
applications where quantile regression is used, understanding the covariate
effect on different quantiles is of paramount interest. Therefore, we seek to
quantify covariate effects on specific aspects of the response distribution as
measured by the quantile function $Q(\tau|\mbox{\bf X})$.

We quantify covariate quantile effects using the accumulative local effects
(ALEs) of Apley and Zhu (2020). The sensitivity of
$Q(\tau|\mathcal{W},\mbox{\bf X})$ to covariate $j$ is naturally quantified by
the partial derivative

$q_{j}(\tau|\mathcal{W},\mbox{\bf X})=\frac{\partial
Q(\tau|\mathcal{W},\mbox{\bf X})}{\partial X_{j}}.$

The ALE begins by averaging $q_{j}(\tau|\mathcal{W},\mbox{\bf X})$ over X
conditioned on $X_{j}=u$,

${\bar{q}}_{j}(\tau|\mathcal{W},u)=\mathbb{E}_{\mbox{\bf
X}}\\{q_{j}(\tau|\mathcal{W},\mbox{\bf X})|X_{j}=u\\}.$

The ALE main effect function of $X_{j}$ is then defined as

$ALE_{j}(\tau|\mathcal{W},x)=\int_{0}^{x}{\bar{q}}_{j}(\tau|\mathcal{W},u)du.$

Analogous formulas define the second-order ALE interaction effect for $X_{j}$
and $X_{l}$, $\mbox{ALE}_{jl}(\tau|\mathcal{W},x_{j},x_{l})$, by taking the
partial derivative with respect to both $X_{j}$ and $X_{l}$. These functions
can be plotted by $\tau$ to summarize how the predicted quantile changes with
respect to change in the covariate values. In addition to the ALEs, we follow
Greenwell et al. (2018) and distill the ALE function to one-number summaries
to compare variable importance by quantile level. The variable importance (VI)
for continuous covariates are characterized by the standard deviation of the
ALE with respect to the marginal distribution of X, i.e.,
$\mbox{VI}_{j}(\tau|\mathcal{W})=\mbox{SD}\\{ALE_{j}(\tau|\mathcal{W},X_{j})\\}$
and
$\mbox{VI}_{jl}(\tau|\mathcal{W})=\mbox{SD}\\{ALE_{jl}(\tau|\mathcal{W},X_{j},X_{l})\\}$.
For discrete covariates with few unique levels, the standard deviation is
replaced by the range.

The ALE and VI summaries depend on the model parameters, $\mathcal{W}$. They
can either be evaluated using a point-estimate $\widehat{\mathcal{W}}$ to give
a point-estimate of the summaries, $ALE_{j}(\tau|\widehat{\mathcal{W}},x)$ and
$\mbox{VI}_{j}(\tau|\widehat{\mathcal{W}})$, or, for a Bayesian analysis, the
posterior samples can be used to quantify uncertainty of the summaries such as
the posterior probability that variable $j$ is more important than variable
$l$ for predicting conditional quantile at $\tau$.

## Computational approaches

SPQR includes four computational algorithms to estimate the model parameters
and importance measures: maximum likelihood estimation (MLE), maximum _a
posteriori_ probability (MAP), Hamiltonian Monte Carlo (HMC) and the no-U-turn
sampler (NUTS). These algorithms are described in detail below.

### Maximum likelihood estimation (MLE)

Directly modeling the conditional PDF means that the negative likelihood
function given $n$ training observations $(\mbox{\bf X}_{i},y_{i})$ for
$i\in\\{1,...,n\\}$ has a closed form,

$\ell(\mathcal{W})=-\sum_{i=1}^{n}\log\left\\{\sum_{k=1}^{K}M_{k}(y_{i})\theta_{k}(\mbox{\bf
X}_{i},\mathcal{W})\right\\}$ (7)

and can be used as a loss function to be minimized, leading to the MLE
estimator

$\widehat{\mathcal{W}}_{\text{MLE}}=\operatorname*{arg\,min}_{\mathcal{W}}\ell(\mathcal{W}).$
(8)

Equation (8) can be solved using standard back-propagation algorithms. The MLE
estimator, however, does not put any constraint on the magnitude of
$\mathcal{W}$ and therefore can result in an unstable model that overfits the
data. One solution is to augment the loss function in (7) with a
regularization term that penalizes large $\mathcal{W}$, for example the weight
decay regularization. However, carefully choosing the penalty coefficients is
crucial to the predictive performance but increases computation complexity.
Furthermore, uncertainty analysis of the estimators obtained by (8) is not
straightforward and has to rely on bootstrap approaches, bringing additional
computational challenges to the already complex problem.

### Bayesian estimation

To address the lack of uncertainty quantification posed by MLE, we adopt a
Bayesian framework and give prior distribution for the weight parameters. The
uncertainty of the SPQR estimators is then characterized by their posterior
distributions. Specifying prior distributions also provides regularization to
the model and stabilizes weight estimation. We assume hierarchical normal
priors for the weights

$W_{ij}^{(l)}|\sigma^{(l)},\lambda_{j}^{(l)}\sim\mathcal{N}\left(0,\sigma^{(l)2}\lambda_{j}^{(l)2}\right),\
p(\lambda_{j}^{(l)})\sim p(\lambda_{j}^{(l)};\gamma_{\lambda})\mbox{\ \ \ for\
\ }j\geq 0$ (9)

where $\sigma^{(l)}$ is the layer-wise global scale shared by all weights in
layer $l$, which can either be set to a constant value or estimated using non-
informative priors, and $\lambda_{j}^{(l)}$ is a unit-wise local scale with
hyper-prior $p(\lambda_{j}^{(l)};\gamma_{\lambda})$. Hierarchical normal
distribution is the most commonly used prior for NN weights as they impose a
weight decay penalty with estimable penalty coefficients, and many priors
proposed in the Bayesian NN literature are variants of (9). In SPQR, we
consider three of such models as summarized in Table 1: the Gaussian Process
(GP) prior, the Automatic Relevance Determinantion (ARD) prior, and the
Gaussian Scale Mixture (GSM) prior.

The GP prior is proposed by Neal (1996) who considers the weights and bias in
each layer as separate parameter blocks. The weights in each layer depend on a
common variance, $\lambda^{(l)}$, and the bias is given a separate variance.
In addition, the variances on weights are scaled by the width of the layer,
$V_{l-1}$, for all but the input layer. The NN under such setting is shown to
converge to a certain multi-output Gaussian process under the condition that
$V_{l}\rightarrow\infty$ for $l\in[1,L-1]$ in both the case of $L=2$ (Neal,
1996) and $L\rightarrow\infty$ (Matthews et al., 2018). The GP prior assumes
that the weights marginally follow Gaussian distributions, since they share a
common variance hyperparameter. Recent studies on Bayesian NNs, however, have
found that distribution of weights in a deep NN (DNN) can be heavy-tailed
(Vladimirova et al., 2019; Fortuin et al., 2021). Therefore, it might be
helpful to incorporate such knowledge and allow a wider model coverage for the
prior distributions for weights. The ARD prior was originally developed by
MacKay (1992) who assigns weights in the input-to-hidden layer unit-wise local
scales so that the magnitude of weights associated with each input will
determine its relevance. Under this setting, the marginal distribution of the
input weights, $W_{ij}^{(1)}|\sigma^{(1)}$, depends on the hyper-prior for
local scales, $\lambda^{(1)}_{j}$, through the integration

$p(w_{ij}^{(1)}|\sigma^{(1)})=\int\mathcal{N}(0,\sigma^{(1)2}\lambda_{j}^{(1)2})p(\lambda_{j}^{(1)};\gamma_{\lambda})d\lambda_{j}^{(1)}.$
(10)

For all other layers, the ARD prior has the same structure as that of the GP
prior. The GSM prior is a direct generalization of the ARD prior by allowing
the layer-wise global scale $\sigma^{(l)}$ to be estimable and all layers to
have the flexibility of (10). Therefore, it not only determines the relevance
of each input feature but also that of each latent feature in deeper layers.

Table 1: Prior distributions: SPQR allows for several models for the prior distribution for the layer-wise global scale, $\sigma^{(l)}$, and unit-wise local scale, $\lambda_{j}^{(l)}$. This table gives the Gaussian Process (GP), Automatic Relevance Determination (ARD) and Gaussian Scale Mixture (GSM) priors in terms of unit index in layer $l$, $j\in[0,V_{l}]$, and hyperparameters $\gamma_{\sigma}$ and $\gamma_{\lambda}$. Name of the prior | $\sigma^{(l)}$ | $\lambda^{(l)}_{0}$ | $\lambda^{(l)}_{j},\ \forall j\geq 1$
---|---|---|---
GP | 1 | $p(\lambda_{0}^{(l)};\gamma_{\lambda})$ | $\lambda^{(l)}\sim p(\lambda^{(l)};\gamma_{\lambda})$
$\lambda_{j}^{(1)}=\lambda^{(1)};\ \lambda_{j}^{(l)}=\lambda^{(l)}/V_{l-1},\
\forall l\geq 2$
ARD | 1 | $p(\lambda_{0}^{(l)};\gamma_{\lambda})$ | $\lambda_{j}^{(1)}\stackrel{{\scriptstyle iid}}{{\sim}}p(\lambda_{j}^{(1)};\gamma_{\lambda})$
$\lambda_{j}^{(l)}\sim\mbox{GP},\ \forall l\geq 2$
GSM | $p(\sigma^{(l)};\gamma_{\sigma})$ | $p(\lambda_{0}^{(l)};\gamma_{\lambda})$ | $\lambda_{j}^{(l)}\stackrel{{\scriptstyle iid}}{{\sim}}p(\lambda_{j}^{(l)};\gamma_{\lambda})$

To complete the prior specification, we assign the local scale,
$\lambda_{j}^{(l)}$, non-informative hyper-priors. In the case of GSM, we also
specify a non-informative hyper-prior for the global scale, $\sigma^{(l)}$. A
common choice is the inverse-Gamma distribution,

$\begin{split}p(\cdot;\gamma_{\lambda})&=\mathcal{IG}amma(a_{\lambda},b_{\lambda})\\\
p(\cdot;\gamma_{\sigma})&=\mathcal{IG}amma(a_{\sigma},b_{\sigma}),\end{split}$
(11)

which is also what will be assumed in this article. Under this setting, the
marginal distribution of weights $p(w_{ij}^{(l)}|\sigma^{(l)})$ follows a
Student-t whose degrees-of-freedom is determined by the hyperparameters,
$a_{\lambda}$ and $b_{\lambda}$. Notice that in the cases of GP prior, by the
scaling property of Gamma distributions, the prior for $\lambda_{j}^{(l)},\
l\geq 2$ is also inverse-Gamma

$\lambda_{j}^{(l)}\stackrel{{\scriptstyle
indep}}{{\sim}}\mathcal{IG}amma(a_{\sigma},\frac{b_{\sigma}}{V_{l-1}})\mbox{\
\ \ for\ \ }j\geq 2$ (12)

#### Maximum a posteriori estimation (MAP)

Before we introduce the fully Bayesian approach to estimate NN parameters
using HMC, we note that the prior specification in (10) enables the adoption
of a MAP method to estimate the weights, i.e.,

$\begin{split}\widehat{\mathcal{W}}_{\text{MAP}},\ \hat{\boldsymbol{\sigma}},\
\hat{\boldsymbol{\Lambda}}&=\operatorname*{arg\,min}_{\mathcal{W},\boldsymbol{\sigma},\boldsymbol{\Lambda}}\left\\{\ell(\mathcal{W})-\log
p(\mathcal{W},\boldsymbol{\sigma},\boldsymbol{\Lambda})\right\\}\\\
p(\mathcal{W},\boldsymbol{\sigma},\boldsymbol{\Lambda})&=\prod_{k=1}^{K}\prod_{l=1}^{L}\prod_{j=0}^{V_{l}}\mathcal{N}(W_{kj}^{(l)}|0,\sigma^{(l)}\lambda_{j}^{(l)})p(\sigma^{(l)}|\gamma_{\sigma})p(\lambda_{j}^{(l)}|\gamma_{\lambda}).\end{split}$
(13)

Solving the optimization problem in (13) will not lead to sensible results as
the gradients of $\sigma^{(l)}$ and $\lambda_{j}^{(l)}$ do not depend on the
data. We adopt a reparameterization of (13) to let the likelihood function
become direct function of the scale hyperparameters

$\begin{split}\widehat{\mathcal{Z}},\ \hat{\boldsymbol{\sigma}},\
\hat{\boldsymbol{\Lambda}}&=\operatorname*{arg\,min}_{\mathcal{Z},\boldsymbol{\sigma},\boldsymbol{\Lambda}}\left\\{\ell(\mathcal{W})-\log
p(\mathcal{Z},\boldsymbol{\sigma},\boldsymbol{\Lambda})\right\\}\\\
W_{kj}^{(l)}&=\sigma^{(l)}\lambda_{j}^{(l)}Z^{(l)}_{kj}\\\
p(\mathcal{Z},\boldsymbol{\sigma},\boldsymbol{\Lambda})&=\prod_{k=1}^{K}\prod_{l=1}^{L}\prod_{j=0}^{V_{l}}\mathcal{N}(Z_{kj}^{(l)}|0,1)p(\sigma^{(l)}|\gamma_{\sigma})p(\lambda_{j}^{(l)}|\gamma_{\lambda}).\end{split}$
(14)

Equation (14) can be solved by substituting the prior distributions for
variance hyperparameters with any of the three models described in Table 1.
The MAP estimates of $\mathcal{W}$ can then be calculated as
$\widehat{W}_{kj,\text{MAP}}^{(l)}=\hat{\sigma}^{(l)}\lambda_{j}^{(l)}\hat{Z}^{(l)}_{kj}$.
As with the MLE method, MAP estimation uses standard back-propagation
algorithms, provides only a point estimate of ${\cal W}$ and does not allow
for direct quantification of model uncertainty.

#### Hamiltonian Monte Carlo (HMC)

MCMC produces samples from the posterior distribution of the parameters that
can be used to approximate their entire posterior distribution. We follow the
strategy of Neal (1996) and use a block-updating scheme that utilizes two MCMC
algorithms. The conditional distribution of weights given the variance
hyperparameters and the data are approximated using the Hamiltonian Monte
Carlo (HMC) sampler, whereas the conditional distribution of variance
hyperparameters given the weights are approximated using the Gibbs sampler
(Geman and Geman, 1984). HMC permits efficient sampling from a high-
dimensional target distribution by using its gradient with respect to each
direction. It uses an approximate Hamiltonian dynamics simulation based on
numerical integration which allows the sampler to explore more carefully in
regions with high density and escape quickly from regions with low density
(Betancourt, 2017). The candidate value of $\mathcal{W}$ found by this
simulation is then accepted with a Metropolis step to correct for any
numerical error resulted from the numerical integration.

We consider two implementations of the HMC sampler: the vanilla HMC sampler
(Neal, 2011) and the more advanced no-U-turn sampling (NUTS; Hoffman and
Gelman, 2014). The vanilla HMC sampler requires setting the number of
approximate integration time $t$, which is the number of leap-frog steps
$L_{\epsilon}$ multiplied by the step size $\epsilon$. The step size is
automatically optimized during warmup sample iterations using dual-averaging
to match an acceptance-rate target, and the number of steps will be calculated
as $L_{\epsilon}=\lfloor t/\epsilon\rfloor$. An Euclidean metric M, where
$\mbox{\bf M}^{-1}$ estimates the posterior covariance of $\mathcal{W}$, is
also estimated during warmup to help project the parameters to a space where
sampling can be done more efficiently. The NUTS, compared to HMC, has the
further advantage of adaptively setting the number of leap-frog steps
$L_{\epsilon}$ on the fly during both warmup and sampling. This greatly
reduces the required effort on the users’ side to select a reasonable value
for the number of approximate integration time $t$.

Let $\mathcal{W}_{s}$ be the posterior samples of $\mathcal{W}$ after
discarding warmup iterations, for $s\in\\{1,...,S\\}$. Since the NN is over-
parameterized and individual weights $W^{(l)}_{ij}\in\mathcal{W}$ are usually
unidentified, the posterior distribution of the weights themselves might not
be very meaningful. However, the samples can produce estimates of meaningful
quantities such as the conditional quantile function (QF) $Q(\tau|\mbox{\bf
X},\mathcal{W})$. For example, the conditional QF estimator is the posterior
mean

$\hat{Q}_{\text{MCMC}}(\tau|\mbox{\bf
X},\mathcal{W})=\frac{1}{S}\sum_{s=1}^{S}Q(\tau|\mbox{\bf X},\mathcal{W}_{s})$

and point-wise credible bands are obtained as the sample quantiles of
$Q(\tau|\mbox{\bf X},\mathcal{W}_{s})$.

Table 2: The overview of functions in package SPQR. Function |  Description
---|---
SPQR() |  Main function of the package. Fits SPQR using the MLE, MAP, or MCMC method. Returns an object of S3 class "SPQR", a list which includes the fitted model (model), the model configuration (config), the control parameters (control), the running time (time), the covariate matrix (X), the response vector (Y), as well as method-dependent training information.
cv.SPQR() |  Fits SPQR using MLE or MAP method, and computes K-fold cross-validation (CV) error.
createFolds.SPQR() |  Generate pre-computed CV folds.
summary() |  Extracts and computes a list of summary information of a "SPQR" class object. Returns an object of S3 class "summary.SPQR".
print.summary() |  Prints the contents of a "summary.SPQR" class object in a user-friendly way.
print() |  Computes and prints the summary information of a "SPQR" class object. Equivalent to print.summary(summary()).
coef() |  Computes and returns the estimated spline coefficients $\theta_{k}(\mbox{\bf X},\widehat{\mathcal{W}})$ of a "SPQR" class object.
predict() |  Computes and returns the estimated PDF/CDF/QF of a "SPQR" class object.
QALE() |  Computes and returns the quantile accumulative local effects (ALE) of a "SPQR" class object.
plotEstimator() |  Computes and plots the estimated PDF/CDF/QF curves of a "SPQR" class object.
plotGOF() |  Performs a visual goodness-of-fit test for the estimated conditional PDF using probability inverse transformation method.
plotMCMCtrace() |  Show trace plot of the log-likelihood or a specified estimate of a "SPQR" class object fitted with method="MCMC".
plotQALE() |  Computes and plots the quantile ALE effects of a "SPQR" class object.
plotQVI() |  Computes and plots the ALE-induced quantile variable importance measures of a "SPQR" class object.
autoplot() |  A wrapper function that creates a user-specified plot for a "SPQR" class object by calling one of the plot functions above.

In summary, the MLE and MAP approaches estimate the parameters by directly
minimizing a loss function and are computationally faster, but they only
provide point estimates of $\mathcal{W}$ as well as the ALE and VI metrics.
HMC is a Markov chain Monte Carlo (MCMC) algorithm that approximates the
posterior distribution of the parameters and is computationally slower, but
produces posterior samples that can be used to quantify uncertainty. However,
HMC requires careful tuning to ensure sampling efficiency. NUTS is an adaptive
variant of HMC that requires little tuning at the further expense of
computational efficiency. Therefore, we recommend the MLE and MAP estimator
for huge datasets where uncertainty quantification is less of a concern, NUTS
for small datasets where uncertainty quantification is critical, and HMC as a
balance between these options.

## The SPQR package

The R package SPQR has two model fitting functions, SPQR() and cv.SPQR(), as
well as various helper functions for handling tasks such as model validation,
model prediction, and results visualization. Table 2 lists all the main
functions provided by SPQR.

### The main fitting function

The SPQR function specifies a semi-parametric conditional density regression
model of the type (5) and estimates the model parameters using one of the four
computational approaches described in the previous section. It has the
following arguments

SPQR(X, Y, n.knots=10, n.hidden=c(10), activation=c("tanh","relu"),
method=c("MLE","MAP","MCMC"), prior=c("GP","ARD","GSM"), hyperpar=list(),
control=list(), normalize=FALSE, verbose=TRUE, seed=NULL, ...)This function
takes two required arguments, a $n\times p$ covariate matrix (without
intercept column) X and a response vector Y. The covariate matrix is expected
to contain only numeric features, and that all categorical features are
converted to numeric values in advance. The covariates are also recommended,
although not required, to be normalized/standardized to have the same scale to
stabilize gradient based optimization that will be used to estimate the
parameters. The response vector, on the other hand, is required to take values
between 0 and 1. We provide a normalize argument such that, when setting
normalize=TRUE, all variables will be scaled to the unit interval using min-
max transformation. Their original scales will be recorded to back-transform
the estimated density and quantile function. By default, however, we set
normalize=FALSE as we want the users to have full control on how variables are
scaled, such as using prior information on the domains of the variables, etc.

The arguments n.knots and n.hidden are the number of basis functions, $K$, and
the number of hidden neurons, $V_{l}$, that define the SPQR model in (5). SPQR
uses the function mSpline() in the splines2 package (Wang and Yan, 2021) to
construct the basis functions. We require setting n.knots to at least 5 as the
model may severely underfit otherwise. The n.hidden argument accepts a vector
of integers such that n.hidden[l] is the the number of hidden neurons in layer
$l$ for $l\in\\{1,...,L-1\\}$. It should be noted that although n.knots and
n.hidden are given default values for convenience, they should be tuned per
application for optimal performance of the SPQR estimators. The argument
activation corresponds to the hidden layer activation function $\phi$ in (3),
and is set to activation="tanh" for hyperbolic tangent by default.

The argument method determines the computational approach to be used to
estimate the model parameters $\mathcal{W}$. The MLE and MAP estimators are
obtained by solving (8) and (13) respectively using gradient-based stochastic
optimization. Specifically, we use the Adam optimizer in the torch package.
The torch package is an R implementation of the open source machine learning
platform – PyTorch (Paszke et al., 2019). It supports hardware acceleration
for systems with a CUDA-compatible NVIDIA graphical processing unit (GPU),
which SPQR takes advantage of. The HMC and NUTS algorithms are implemented
using C++ and mirror those in Stan (Stan Development Team, 2022). We did not
use Stan directly for two reasons. First, Stan uses automatic differentiation
whereas the analytical gradients of NNs are fairly easy to derive and
evaluate. SPQR depends on Rcpp and RcppArmadillo to efficiently compute the
log-posterior and its gradients. Secondly, Stan does not allow block-updating
using both HMC and Gibbs sampler that exploits the conjugacy of hyper-priors.
By default, we set method="MAP" since it is computationally faster than the
MCMC method and less prone to overfit than the MLE method.

The argument prior is used only for the Bayesian methods and corresponds to
one of the three variance hyperpriors described in Table 1. We set prior="GP"
as default for simplicity. The argument hyperpar is a list of named hyper-
prior hyperparameters to use instead of the default values, including
a_lambda, b_lambda, a_sigma and b_sigma. They correspond to $a_{\lambda}$,
$b_{\lambda}$, $a_{\sigma}$ and $b_{\sigma}$ in (11). The default value is
0.001 for all four hyperparameters.

The argument control is a list of named and method-dependent parameters that
allows finer control of the behavior of the computational approaches. The
available parameters for MLE and MAP estimators are shown in Table 3. The NNs
for MLE and MAP estimators are structured using a templated module that
consists of only fully-connected layers ("nn_linear"), batch normalization
("nn_batch_norm1d"), and dropout ("nn_dropout"). Dropout and batch
normalization are not used by default, but may be useful when training deep
and wide NNs. We recommend setting use.GPU=TRUE to further accelerate
computation on CUDA-configured systems. Early stopping is implemented to avoid
overfitting. We use valid.pct$\times$100% of the data as the validation set.
During each epoch, a snapshot of the trained model is saved in
save.path/save.name if it leads to a decrease in the validation loss. When the
validation loss stops decreasing for early.stopping.epochs, the training stops
and the best model is loaded and returned.

Table 3: Control parameters for MLE and MAP. Parameter |  Description
---|---
use.GPU |  A Boolean flag for GPU utilization. Default is FALSE.
lr |  Learning rate of Adam optimizer. Default is 0.01.
dropout |  Dropout probabilities. A length-two vector of which the first entry specifies the dropout probability in the input-to-hidden layer and the second entry specifies the probabilities in all hidden-to-hidden layers. Default is c(0,0) which corresponds to no dropout.
batchnorm |  A Boolean flag for batch normalization. Default is FALSE.
epochs |  The (maximum) number of passes of the entire training dataset by Adam. Default is 200.
batch.size |  Size of mini batches to calculate gradient. Default is 128.
valid.pct |  Percentage of data used as validation set. Default is 0.2.
early.stopping.epochs |  The number of epochs before stopping if the validation loss does not decrease. Default is 10.
print.every.epochs |  The number of epochs before next training progress in printed. Default is 10.
save.path |  Path to save the trained torch model with the lowest validation loss. Default is file.path(getwd(),"SPQR_model").
save.name |  File name to save the trained torch model with the lowest validation loss. Default is "SPQR.model.pt".

The available parameters for MCMC estimator are shown in Table 4. These
parameters have similar meanings to those in the stan() function in the rstan
package (Guo et al., 2021), and detailed explanations of their effects can be
found in the Stan reference manual (Stan Development Team, 2022). By default,
SPQR uses algorithm="NUTS" to approximate the posterior distribution of
$\mathcal{W}$ as it adaptively selects both the leap-frog discretization step-
size $\epsilon$ as well as the number of steps $L_{\epsilon}$ which are
crucial to the sampling efficiency of HMC. The value of step-size is also
affected by the target Metropolis acceptance rate delta. When the NN is large
and its posterior geometry is complex, a larger delta allows NUTS to explore
the posterior more carefully using smaller steps. However, a smaller step-size
may require a larger max.treedepth to avoid premature stopping which will
significantly increase the algorithm run-time.

Finally, the argument verbose determines whether training progress should be
printed, seed sets the seed for random number generation when reproducibility
is desired, and ... allows any of the control parameters to be specified
directly in the function call instead of in control.

Table 4: Control parameters for MCMC. These parameters are similar to those in stan() in rstan. Detailed explanations can be found in the Stan reference manual.  Parameter |  Description
---|---
algorithm |  The sampling algorithm; "HMC": Hamiltonian Monte Carlo with dual-averaging, "NUTS": No-U-Turn sampler (default).
iter |  The number of iterations (including warmup). Default is 2000.
warmup |  The number of warmup/burn-in iterations for step-size and mass matrix adaptation. Default is 500.
thin |  The period for saving post-warmup samples. The default is 1.
stepsize |  The discretization interval/step-size $\epsilon$ of leap-frog integrator. Default is NULL which means it will be adaptively selected during warmup iterations.
metric |  The mass matrix M; "unit": diagonal matrix of ones, "diag": diagonal matrix with positive diagonal entries estimated during warmup iterations (default), "dense": a dense, symmetric positive definite matrix during warmup iterations.
delta |  The target Metropolis acceptance rate. Default is 0.9.
max.treedepth |  The maximum tree depth in NUTS. The default is 6.
int.time |  The integration time $t$ in HMC. The number of leap-frog steps is then calculated as $L_{\epsilon}=\lfloor t/\epsilon\rfloor$. Default is 1.

### Cross-validation function

As seen above, the Adam routine used to compute MLE and MAP estimators has
numerous potential tuning parameters, such as the learning rate and batch
size. In addition, the number of basis functions and hidden neurons will also
affect the quality of the estimator. To allow users to select the best model
by comparing different model configurations, we provide the cv.SPQR() function
that calculates the cross-validation (CV) error for a given configuration of
the MLE or MAP estimator. In addition to all arguments of SPQR, cv.SPQR()
takes the argument folds which are pre-computed folds on which CV can be
performed. It should be noted that cv.SPQR() itself does not perform model
selection, but can be used as a building block in a grid search loop. There
are many ways to generate pre-computed CV folds, such as the createFolds()
function in the caret package (Kuhn, 2022). For users’ convenience, we provide
the function createFolds.SPQR() which is equivalent to createFolds() function
in caret but returns an unnamed list. Among all control parameters listed in
Table 3, we recommend tuning lr, in addition to n.knots and n.hidden. A
reasonable range of values to consider is n.knots $\in\\{8,10,12\\}$, n.hidden
$\in\\{8,10,15\\}$ and lr $\in\\{e^{-6},e^{-5},e^{-4},e^{-3}\\}$.

The cv.SPQR() function is only applicable to MLE and MAP estimators. For the
MCMC estimator, prediction accuracy can be measured by the expected log
pointwise predictive density (ELPD) which can be estimated by Bayesian leave-
one-out (LOO) CV or by widely applicable information criterion (WAIC, Vehtari
et al., 2017). Calculation of these statistics is integrated in the summary()
functionality and described in the following section. Once convergence is
achieved, we recommend running separate chains with different combinations of
n.knots and n.hidden and select the best model with the highest ELPD.

### Helper functions for "SPQR" object

The SPQR() function returns an object of S3 class "SPQR", a compound list that
contains the fitted model besides various other information. The summary()
function summarizes the output produced by SPQR() and structures them in a
more organized way to be examined by the user. For SPQR() fitted using any of
the four methods, the output returned by summary() contains the estimation
method, the run-time, and the NN architecture. For MLE and MAP estimators, the
output contains the training and validation loss of the final (best if early
stopping is used) model, and selected information about the optimizer such as
learning rate and batch size. For MAP and MCMC estimators, the output contains
the variance hyper-prior model. For MCMC estimator, summary() calculates
various diagnostic statistics that can be used to evaluate the fit of the
model and convergence of the MCMC chain. Specifically, we calculate Bayesian
LOO-CV and WAIC using the loo package (Vehtari et al., 2020) for model
comparison, and the average Metropolis acceptance ratio and the number of
divergences of post-warmup iterations to determine if the chain is reliable.
The output of summary() is also an object of S3 class "summary.SPQR", to which
print.summary() can be applied to print the aforementioned contents in a user-
friendly way. The function print.summary() has an optional argument showModel
that when set to TRUE, additionally prints the NN architecture by layer.
Instead of using summary and print.summary, the user can achieve the same goal
by directly print() the "SPQR" object in one function call.

The coef() function outputs the estimated spline coefficients
$\theta_{k}(\mbox{\bf X},\widehat{\mathcal{W}})$ given X. For MCMC estimator,
these will be the posterior means of $\theta_{k}(\mbox{\bf X},\mathcal{W})$.
The predict() function computes different estimates of the conditional
distribution. When type="PDF" or type="CDF", it computes $f(Y|\mbox{\bf
X},\widehat{\mathcal{W}})$ or $F(Y|\mbox{\bf X},\widehat{\mathcal{W}})$ for
every combination of X and Y; when type="QF", it computes $Q(\tau|\mbox{\bf
X},\hat{\mathcal{W}})$ for every combination of X and tau. The argument Y is
optional, and when left unspecified predict() will use nY to define a grid on
[0,1] for estimation. This is useful when the user wants to estimate the full
PDF/CDF/QF curve. For MCMC estimator, two additional arguments are available:
ci.level and getAll. The argument ci.level allows the user to extract
ci.level$\times$100% credible intervals of the estimates in addition to the
posterior means, whereas setting getAll=TRUE allows the user to extract all
posterior samples of the corresponding estimates.

### Quantile accumulated local effect (QALE)

The function QALE() is largely based on the ALEPlot() function in the ALEPlot
package (Apley, 2018), but adapted for the quantile regression setting to
compute ALEs at different quantiles. The function takes the following
arguments,

QALE(object, var.index, tau=seq(0.1,0.9,0.1), n.bins=40, ci.level=0,
getAll=FALSE, pred.fun=NULL).The argument object corresponds to the fitted
SPQR object of class "SPQR". The var.index is a numeric scalar or length-two
vector of indices of the covariates for which the ALEs will be calculated.
When length(var.index) = 1, the function calculates the main effects for
X[,var.index]; when length(var.index) = 2, the function calculates interaction
effects for X[,var.index[1]] and X[,var.index[2]]. The argument tau is a
numeric vector of quantile levels at which the ALEs will be calculated. The
argument n.bins is a numeric scalar that specifies the maximum number of
intervals into which the covariate range is divided when calculating the ALEs.
The actual number of intervals depends on the number of unique values in
X[,var.index]. When length(var.index) = 2, n.bins is applied to both
covariates. The arguments ci.level and getAll allow uncertainty analysis of
the calculated ALEs, but are only implemented for main effects for a single
covariate. Finally, the argument pred.fun accepts a custom quantile prediction
function that will be used instead of the built-in predict() function for
calculating ALE. This can be useful when the user wants to compare the QALE
calculated using SPQR to that using other quantile regression models, or maybe
that using the true model in a simulation study.

### Plot functions

Various plot functions are implemented to visualize results from model
prediction, model diagnostics and model interpretation. All of them take
directly the "SPQR" object as input and return a "ggplot" object that can be
further customized using layers from the ggplot2 package (Wickham, 2016). The
plotEstimator() function computes and plots the estimated PDF/CDF/QF curve for
a single observation. The plotGOF() function plots the quantiles of
probability integral transform of the observed data against the quantiles of a
uniform distribution. Let $U_{i}=F(Y_{i}|\mbox{\bf
X}_{i},\widehat{\mathcal{W}})$ be the estimated CDF of the $i$th observation
in the dataset. By the probability integral transform, if the observed data
indeed distribute according to $F(Y|\mbox{\bf X},\widehat{\mathcal{W}})$ then
the sample $U_{i},\ i\in\\{1,...,n\\}$ should correspond to independent
samples of a standard uniform distribution, i.e.,
$U_{1},...,U_{n}\stackrel{{\scriptstyle iid}}{{\sim}}\mathcal{U}(0,1)$. Thus
the Quantile-Quantile plot (Q-Q plot) created by plotGOF() can be used for
visually checking the goodness-of-fit of the SPQR estimator. When SPQR is
fitted with method = "MCMC", the function plotMCMCtrace() can be used to show
trace plot of a target to examine the convergence and autocorrelation of the
MCMC chain. Available target are "loglik" for the log-likelihood, and "PDF",
"CDF" and "QF" for corresponding estimate for a single observation. The
function also takes an optional argument window which allows examination of
different parts of the chain. The functions plotQALE() and plotQVI() support
the QALE() function in helping the user understand how the covariates affect
different quantiles, and thus accept similar arguments. In particular, the
function plotQALE() plots the estimated quantile ALE main effects when
length(var.index) = 1 and quantile ALE interaction effects when
length(var.index) = 2, whereas the function plotQVI() compares the quantile
ALE-induced variable importance (VI) between all considered covariates. The
var.index argument can be left unspecified in plotQVI() in which case all
covariates in the data set are considered. Similar to the predict() function,
most plot functions accept arguments ci.level or getAll (or both) for
uncertainty quantification when SPQR is fitted with method = "MCMC". Table 5
summarizes the implementation of uncertainty quantification for plot functions
in SPQR. Finally, the function autoplot() provides a wrapper for all the plot
functions mentioned above.

Table 5: Implementation of uncertainty quantification for plot functions in SPQR. Name of the function | ci.level | getAll
---|---|---
plotEstimator() | $\checkmark$ | $\checkmark$
plotGOF() | $\times$ | $\checkmark$
plotMCMCtrace() | N/A | N/A
plotQALE() | $\checkmark$ | $\checkmark$
plotQVI() | $\checkmark$ | $\times$

## Examples

### Simulation

We start with a simple simulation study to demonstrate the effectiveness of
SPQR in estimating conditional density and quantile functions, as well as
usages of functions introduced in the previous section. We consider a three
dimensional covariate $\mbox{\bf X}=(X_{1},X_{2},X_{3})$ with variables
independently generated from a uniform distribution. The response $Y$ follows
a Beta distribution whose shape parameters are functions of $X_{1}$ and
$X_{2}$, the third covariate is irrelevant.

$\begin{split}X_{j}&\stackrel{{\scriptstyle iid}}{{\sim}}\mathcal{U}(0,1),\
j=1,...,3\\\
Y&\sim\mathcal{B}eta\left(10\mbox{expit}\left\\{1-5X_{1}X_{2}\right\\},10\left[1-\mbox{expit}\left\\{1-5X_{1}X_{2}\right\\}\right]\right).\end{split}$
(15)

Here $\mbox{expit}(u)=1/(1+e^{-u})$ is the inverse logistic link function.
Based on the definition of Beta distribution, the following properties of the
conditional distribution $f(Y|\mbox{\bf X})$ can be obtained

$\begin{split}\mathbb{E}(Y|\mbox{\bf
X})&\propto\mbox{expit}\left\\{1-5X_{1}X_{2}\right\\}\\\
\mbox{Var}(Y|\mbox{\bf
X})&\propto\mbox{expit}\left\\{1-5X_{1}X_{2}\right\\}\left[1-\mbox{expit}\left\\{1-5X_{1}X_{2}\right\\}\right]\\\
\mbox{Skewness}(Y|\mbox{\bf
X})&\propto-\frac{\sqrt{\mbox{expit}\left\\{1-5X_{1}X_{2}\right\\}}}{\sqrt{1-\mbox{expit}\left\\{1-5X_{1}X_{2}\right\\}}}.\end{split}$
(16)

That is, the location, scale and shape of $f(Y|\mbox{\bf X})$ all depend on X.
Figure 1 shows the true density and quantile functios of $Y$ as conditioned on
different combinations of $X_{1}$ and $X_{2}$. We can clearly observe a
varying effect of X on the density and quantile functions of $Y$.

Figure 1: True conditional density and quantile functions of (15) for
different combinations of $X_{1}$ and $X_{2}$.

We generate $n=1000$ samples from (15) and fit SPQR using the MLE, MAP, and
MCMC methods. For all three methods, the default configuration of n.knots=10,
n.hidden=c(10) and activation="tanh" are used. For the MLE and MAP methods,
the learning rate is selected from $\\{e^{-6},e^{-5},e^{-4},e^{-3}\\}$ using
5-fold cross-validation, and the gradients are calculated using mini-batches
of size 256 with a maximum training time of 500 epochs. For the MAP method,
the default prior="ARD" is used.

## Generate simulated data set.seed(919) n <\- 1000 p <\- 3 X <\-
matrix(runif(n*p),n,p) expit <\- 1/(1+exp(-1+5*X[,1]*X[,2])) Y <\-
rbeta(n,10*expit,10*(1-expit)) ## Specifying MLE control parameters
mle.control <\- list(batch.size = 256, epochs = 500, use.GPU = TRUE) ##
Generate CV folds folds <\- createFolds.SPQR(Y, nfold=5) ## Select optimal
learning rate using 5-fold CV lr.grid <\- exp(-6:-3) cve <\- sapply(lr.grid,
FUN=function(lr) { mle.control$lr <\- lr cv.out <\- cv.SPQR(X=X, Y=Y,
folds=folds, method="MLE", control=mle.control) cv.out$cve }) mle.control$lr
<\- lr.grid[which.min(cve)] mle.fit <\- SPQR(X=X, Y=Y, method="MLE",
control=mle.control) ## Specifying MAP control parameters map.params <\-
mle.params cve <\- sapply(lr.grid, FUN=function(lr) { map.control$lr <\- lr
cv.out <\- cv.SPQR(X=X, Y=Y, folds=folds, method="MAP", control=map.control)
cv.out$cve }) map.control$lr <\- lr.grid[which.min(cve)] map.fit <\- SPQR(X=X,
Y=Y, method="MAP", control=map.control)For the MCMC method, the posterior
distribution is approximated using NUTS. We run NUTS for a total of 1000
iterations, discard the first 250 as warm-ups, and save every other iteration.

## Specifying MCMC control parameters mcmc.control <\- list(iter = 1000,
warmup = 250, thin = 2) mcmc.fit <\- SPQR(X=X, Y=Y, method="MCMC",
control=mcmc.control)The function print() returns a short summary of the
results of the fitted object. The summary depends on the fitted method. Here
we show summary for models fitted with method="MLE" and method="MCMC". The
summary for method="MAP" is mostly the same as that for method="MLE" and thus
omitted.

mle.fit#> SPQR fitted using MLE approach#>#> Learning rate: 0.04978707#> Batch
size: 256#>#> Loss:#> train = -140.2645, validation = -142.1546#>#> Elapsed
time: 0.22 minutesmcmc.fit#> SPQR fitted using MCMC approach with ARD
prior#>#> MCMC diagnostics:#> Final acceptance ratio is 0.91 and target is
0.9#>#> Expected log pointwise predictive density (elpd) estimates:#> elpd.LOO
= 683.4493, elpd.WAIC = 683.9026#>#> Elapsed time: 3.05 minutesAn overview of
the NN architecture can be additionally printed by setting the argument
showModel to TRUE.

print(mle.fit, showModel = TRUE)#> SPQR fitted using MLE approach#>#> Learning
rate: 0.04978707#> Batch size: 256#>#> Model specification:#> Layers#> Input
Output Activation#> 3 10 tanh#> 10 10 softmax#>#> Loss:#> train = -140.2645,
validation = -142.1546#>#> Elapsed time: 0.22 minutesThe function plotGOF()
uses a Q-Q plot to inspect the alignment of the probability integral transform
(PIT) and the uniform distribution, which can be used to visually check the
goodness of fit of SPQR(). Figure 2 shows the goodness of fit of the three
SPQR estimators. For SPQR fitted with method="MCMC", the argument getAll is
used to additionally plot all posterior samples of the Q-Q plot. The plots
show that the distribution of PIT is most similar to the uniform distribution
for the MCMC estimator, suggesting that the model based on the MCMC method
best fits the observed data.

## Q-Q plot of PIT for goodness of fit test plotGOF(mle.fit) + ggtitle("Q-Q
plot, MLE") plotGOF(map.fit) + ggtitle("Q-Q plot, MAP") plotGOF(mcmc.fit,
getAll = TRUE) + ggtitle("Q-Q plot, MCMC")

Figure 2: Examples of plot produced by plotGOF(). Goodness-of-fit test for
SPQR estimators.

The reliability of the MCMC estimator can be further assessed by the function
plotMCMCtrace(), which plots traceplots that can be used to visually inspect
the sampling behavior and convergence of the post-warmup chain. The code below
provides two examples of such usage, one uses target="loglik" to plot the
traceplot of the log-likehood; the other uses target="QF", X=c(0.5,0.5,0.5)
and tau=0.5 to plot the traceplot of the estimated median when $\mbox{\bf
X}=(0.5,0.5,0.5)^{\intercal}$. The results are shown in Figure 3 and provide
visual evidence that the chain converged.

## Traceplot of log-likelihood plotMCMCtrace(mcmc.fit, target = "loglik") ##
Traceplot of estimated quantile plotMCMCtrace(mcmc.fit, target = "QF", X =
c(0.5,0.5,0.5), tau = 0.5)

Figure 3: Examples of traceplot produced by plotMCMCtrace().

To investigate the numerical performance of the fitted SPQR models, we use the
function predict() to compute the predicted conditional PDFs and QFs for 3
out-of-sample observations and compare them with their true values.

## Prediction based on 3 out-of-sample observations set.seed(3) X_test <\-
matrix(runif(3*3), nrow = 3) ## PDFs ## Not specifying `Y` means to estimate
on a fine grid ## Default is Y <\- seq(0,1,length=101) pdf.mle <\-
predict(mle.fit, X = X_test, type = "PDF") pdf.map <\- predict(map.fit, X =
X_test, type = "PDF") pdf.mcmc <\- predict(mcmc.fit, X = X_test, type = "PDF")
## QFs qf.mle <\- predict(mle.fit, X = X_test, type = "QF") qf.map <\-
predict(map.fit, X = X_test, type = "QF") qf.mcmc <\- predict(mcmc.fit, X =
X_test, type = "QF") ## Function to compute true PDFs dYgivenX <\- function(y,
X){ expit <\- 1/(1+exp(-1+5*X[1]*X[2])) out <\- dbeta(y, 10*expit,
10*(1-expit)) return(out)} ## Compare the predicted PDFs with their true
values yyy <\- seq(0, 1, length.out = 101) par(mfrow = c(2,3)) for(i in 1:3){
pdf0 <\- dYgivenX(yyy, X_test[i,]) # True pdf1 <\- pdf.mle[i,] # MLE pdf2 <\-
pdf.map[i,] # MAP pdf3 <\- pdf.mcmc[i,] # MCMC plot(yyy, pdf0,
ylim=c(0,1.2*max(pdf0)), type="l", xlab="Y", ylab="PDF",
main=paste("Observation",i), cex.lab=1.2, cex.axis=1.2) lines(yyy, pdf1,
col=2) lines(yyy, pdf2, col=3) lines(yyy, pdf3, col=4) if(i == 1) {
legend("topleft", c("True","MLE","MAP","MCMC"), lty=1, col=1:4, bty="n") }} ##
Function to compute true QFs qYgivenX <\- function(tau,X){ expit <\-
1/(1+exp(-1+5*X[1]*X[2])) out <\- qbeta(tau, 10*expit, 10*(1-expit))
return(out)} ## Compare the predicted QFs with their true values for(i in
1:3){ qf0 <\- qYgivenX(tau,X_test[i,]) # True qf1 <\- qf.mle[i,] # MLE qf2 <\-
qf.map[i,] # MAP qf3 <\- qf.mcmc[i,] # MCMC
plot(tau,qf0,ylim=c(0,1.2*max(qf0)),type="l",
xlab="tau",ylab="Quantile",main=paste("Observation",i), cex.lab=1.2,
cex.axis=1.2) lines(tau,qf1,col=2) lines(tau,qf2,col=3) lines(tau,qf3,col=4)
}The estimated and true PDFs are shown in the top row of Figure 4, and the
estimated and true QFs are shown in the bottom row.

Figure 4: Estimated and true PDFs (top) and QFs (bottom) for 3 out-of-sample
observations.

For the MCMC estimator, we also compute the 95% pointwise credible intervals
using the argument ci.level=0.95 and plot them in Figure 5. These plots show
that SPQR successfully captured the varying effect of X on the distribution of
$Y$.

## PDF point estimates and credible intervals pdf.mcmc <\- predict(mcmc.fit,
X_test, type="PDF", ci.level=0.95) ## QF point estimates and credible
intervals qf.mcmc <\- predict(mcmc.fit, X_test, type="QF", ci.level=0.95) ##
Compare the predicted PDFs with their true values par(mfrow=c(2,3)) for (i in
1:3) { pdf0 <\- dYgivenX(yyy,X_test[i,]) # True pdf1.lb <\-
pdf.mcmc[i,,"lower.bound"] pdf1 <\- pdf.mcmc[i,,"mean"] pdf1.ub <\-
pdf.mcmc[i,,"upper.bound"] plot(yyy,pdf0,ylim=c(0,1.5*max(pdf0)),type="l",
xlab="Y",ylab="PDF",main=paste("Observation",i), cex.lab=1.2, cex.axis=1.2)
lines(yyy,pdf1,col=2) lines(yyy,pdf1.lb,col=2,lty=2)
lines(yyy,pdf1.ub,col=2,lty=2) if (i==1) {
legend("topright",c("True","MCMC"),lty=1:2,col=1:2,bty="n") } } ## Compare the
predicted QFs with their true values for (i in 1:3) { qf0 <\-
qYgivenX(tau,X_test[i,]) # True qf1.lb <\- qf.mcmc[i,,"lower.bound"] qf1 <\-
qf.mcmc[i,,"mean"] qf1.ub <\- qf.mcmc[i,,"upper.bound"]
plot(tau,qf0,ylim=c(0,1.5*max(qf0)),type="l",
xlab=expression(tau),ylab="Quantile",main=paste("Observation",i), cex.lab=1.2,
cex.axis=1.2) lines(tau,qf1,col=2) lines(tau,qf1.lb,col=2,lty=2)
lines(tau,qf1.ub,col=2,lty=2) }

Figure 5: Estimated and true PDFs (top) and QFs (bottom), along with 95%
credible bands, for 3 out-of-sample observations.

Figures 4 and 5 can be easily reproduced using plotEstimator() which inherits
most arguments of predict(). The codes below show examples of plotting the
predicted PDF along with 95% credible bands and plotting the predicted QF
along with all posterior samples. The results are shown in Figure 6.

## PDF with 95% credible bands plotEstimator(mcmc.fit, X = X_test[2,], type =
"PDF", ci.level = 0.95) ## QF with all posterior samples
plotEstimator(mcmc.fit, X = X_test[2,], type = "QF", getAll = TRUE)

Figure 6: Examples of plot produced by plotEstimator(). Left: estimated PDF
along with 95% credible bands. Right: estimated QF along with all posterior
samples.

In most applications where quantile regression is used, understanding the
covariate effect on the conditonal quantile/distribution is of paramount
importance. The function QALE() can be used to quantify either the main effect
of a single covariate or the interaction effect between two covariates on the
predicted quantiles. the codes below show examples of computing the QALE for
$X_{1}$, $X_{2}$, $X_{3}$ respectively at $\tau=0.25$ using fitted "SPQR"
objects, as well as an example using custom prediction function (in this case
it is the true quantile function). To compute the QALE using a custom
prediction function, the user should first define a function that takes X, the
covariate matrix, and tau, a vector of quantile levels, as inputs and returns
an nrow(X) by length(tau) matrix of predicted quantiles. The user should then
pass the defined function to the pred.fun argument and list(X=X) to the object
argument. The results are shown in Figure 7 and provide visual evidence that
SPQR accurately captures the quantile covariate effects.

## Quantile of interest tau <\- 0.25 ## A custom prediction function that
computes the true quantiles pred.fun <\- function(X, tau) { out <\-
matrix(nrow=nrow(X), ncol=length(tau)) for (i in 1:nrow(X)) { out[i,] <\-
qYgivenX(tau, X[i,]) } return(out) } ## Main effect for X1, X2, X3
par(mfrow=c(1,3)) for (j in 1:3) { ale.mle <\- QALE(mle.fit, var.index=j,
tau=tau) # MLE ale.map <\- QALE(map.fit, var.index=j, tau=tau) # MAP ale.mcmc
<\- QALE(mcmc.fit, var.index=j, tau=tau) # MCMC ## The following line will
work as long as `object` is a list that contains `X` ale.ans <\-
QALE(list(X=X), var.index=j, tau=tau, pred.fun=pred.fun) # True
plot(ale.ans$x, ale.ans$ALE, type="l", xlab=parse(text=paste0("X[",j,"]")),
ylab="ALE", cex.lab=1.2, cex.axis=1.2) lines(ale.mle$x, ale.mle$ALE, col=2)
lines(ale.map$x, ale.map$ALE, col=3) lines(ale.mcmc$x, ale.mcmc$ALE, col=4) if
(j == 1) { legend("topright", c("True","MLE","MAP","MCMC"), lty=1, col=1:4,
bty="n") } }

Figure 7: Quantile accumulative local effects (ALEs) for $X_{1}$, $X_{2}$ and
$X_{3}$ respectively at $\tau=0.5$.

We also compute the QALE interaction effect between $X_{1}$ and $X_{2}$ at
$\tau=0.5$. The results are shown in Figure 8. The plots show that SPQR is
able to recover the complex interaction effect between $X_{1}$ and $X_{2}$.

par(mfrow=c(2,2)) ## Interaction effect between X1 and X2 ale.mle <\-
QALE(mle.fit, var.index=c(1,2), tau=0.5) ale.map <\- QALE(map.fit,
var.index=c(1,2), tau=0.5) ale.mcmc <\- QALE(mcmc.fit, var.index=c(1,2),
tau=0.5) ale.ans <\- QALE(list(X=X), var.index=c(1,2), tau=0.5,
pred.fun=pred.fun) ## Contour plots for visualizing 2D effects
image(ale.ans$x[[1]], ale.ans$x[[2]], ale.ans$ALE[,,1], xlab =
parse(text="X[1]"), ylab = parse(text="X[2]"), main = "True", cex.lab=1.1,
cex.axis=1.1) image(ale.mle$x[[1]], ale.mle$x[[2]], ale.mle$ALE[,,1], xlab =
parse(text="X[1]"), ylab = parse(text="X[2]"), main = "MLE", cex.lab=1.1,
cex.axis=1.1) image(ale.map$x[[1]], ale.map$x[[2]], ale.map$ALE[,,1], xlab =
parse(text="X[1]"), ylab = parse(text="X[2]"), main = "MAP", cex.lab=1.1,
cex.axis=1.1) image(ale.mcmc$x[[1]], ale.mcmc$x[[2]], ale.mcmc$ALE[,,1], xlab
= parse(text="X[1]"), ylab = parse(text="X[2]"), main = "MCMC", cex.lab=1.1,
cex.axis=1.1)

Figure 8: Quantile accumulative local effects (ALEs) interaction effect
between $X_{1}$ and $X_{2}$ at $\tau=0.5$.

Figures 7 and 8 can be reproduced using the function plotQALE(). the codes
below show two examples. The first one plots the QALE main effect of $X_{1}$
at $\tau\in\\{0.1,0.5,0.9\\}$, along with 95% credible bands. The second one
plots the QALE interaction effect between $X_{1}$ and $X_{2}$ at
$\tau\in\\{0.1,0.5,0.9\\}$. The results are shown in Figures 9 and 10,
respectively. The results suggest that the main effect of $X_{1}$ and the
interaction effect between $X_{1}$ and $X_{2}$ are more prominent on upper
quantiles than on lower quantiles.

## Main effect for X1 plotQALE(mcmc.fit, var.index = 1, tau = c(0.1,0.5,0.9),
ci.level = 0.95) ## Interaction effect between X1 and X2 plotQALE(mcmc.fit,
var.index = c(1,2), tau = c(0.1,0.5,0.9))

Figure 9: Example of plot produced by plotQALE(). QALE main effect for $X_{1}$
at $\tau\in\\{0.1,0.5,0.9\\}$. Figure 10: Example of plot produced by
plotQALE(). QALE interaction effect between $X_{1}$ and $X_{2}$ at
$\tau\in\\{0.1,0.5,0.9\\}$.

The function plotQVI() computes the QALE-induced variable importance (VI)
measures of each covariate. The covariates are then ranked accordingly in a
barplot at each quantile of interest. the codes below show an example of
computing the quantile VI for $X_{1}$, $X_{2}$ amd $X_{3}$ at
$\tau\in\\{0.1,0.5,0.9\\}$, respectively. The argument ci.level = 0.95 is used
to plot 95% error bar for each VI measure. The result is shown in Figure 11.
The result suggests that $X_{1}$ and $X_{2}$ have similar and significant
effect on quantiles of $Y$ whereas $X_{3}$ has no effect. The effects of
$X_{1}$ and $X_{2}$ are also more prominent on upper quantiles than on lower
quantiles. These observations match the truth.

## Quantile variable importance with error bars plotQVI(mcmc.fit,
tau=c(0.1,0.5,0.9), ci.level=0.95)

Figure 11: Example of plot produced by plotQVI(). Quantile variable importance
(VI) at $\tau\in\\{0.1,0.5,0.9\\}$.

### Australia electricity demand data

For real data application we consider the electricity demand data from Sydney,
Australia. The dataset contains electric energy consumption, recorded from
smart meters, of 247 anonymized residential customers during the period
between July 3rd 2010 and June 30th 2011. It is available from the R package
qgam and has been analyzed by Fasiolo et al. (2021) in the context of quantile
regression. The data can be loaded using the following codes,

## Load the electricity demand data data("AUDem", package = "qgam") meanDem
<\- AUDem$meanDem head(meanDem) doy tod dem dow temp date dem48 #> 1 184 18.0
0.8248777 Sat 3.357407 2010-07-03 17:30:00 0.8978636 #> 2 184 18.5 0.8686110
Sat 2.517073 2010-07-03 18:00:00 0.9417633 #> 3 184 19.0 0.8519471 Sat
1.898399 2010-07-03 18:30:00 0.9148921 #> 4 184 19.5 0.8457030 Sat 1.666667
2010-07-03 19:00:00 0.8982188 #> 5 184 20.0 0.8464793 Sat 1.943989 2010-07-03
19:30:00 0.8942396 #> 6 184 20.5 0.8471040 Sat 2.679803 2010-07-03 20:00:00
0.8818705and the definitions of the variables are given as follows

  * •

doy: the day of the year, from 1 to 365;

  * •

tod: the time of day, ranging from 18 to 22, where 18 indicates the period
from 17:00 to 17:30, 18.5 the period from 17:30 to 18:00 and so on;

  * •

dem: the demand (in kWh) during a 30min period, averaged over the 247
households;

  * •

dow: factor variable indicating the day of the week;

  * •

temp: the external temperature at Sydney airport, in degrees Celsius;

  * •

date: local date and time;

  * •

dem48: the lagged mean demand, that is the average demand (dem) during the
same 30min period of the previous day.

In particular, we are interested in analyzing the effect of temperature and
time on the average demand distribution and its different quantiles. Hence the
response variable is chosen to be dem, and the set of covariates contains doy,
tod, temp and dummy variables representing dow.

library(tidyr) ## Encode dow using dummy variables meanDem <\- meanDem %>%
mutate(value=1) %>% spread(dow, value, fill=0) Y <\- meanDem$dem X <\-
as.matrix(meanDem %>% select(-c(date,dem,dem48)))Given that the sample size is
fairly large, we consider fitting SPQR using the MLE method. The modeling
parameters are selected using a grid search. Specifically, we select the
number of basis functions (n.knots) from $\\{10,15,20\\}$ and the number of
hidden neurons (n.hidden) from $\\{10,15,20\\}$. We focus on 2-hidden-layer
neural networks so n.hidden represents the number of neurons in each hidden
layer. For the control hyperparameters, we also select learning rate (lr) from
$\\{e^{-3},e^{-4},e^{-5}\\}$; we set the batch size (batchsize) to be 128,
maximum number of epochs (epochs) to be 800, and early stopping criterion
(early.stopping.epochs) to be 50. The best model configuration is selected
based on 10-fold CV error.

To accelerate the grid search routine, the function foreach() from the foreach
package (Microsoft and Weston, 2022) can be used to run jobs in parallel.
However, since tensor computation in the torch package is very memory-
consuming, using many cores can easily exceeds the memory limit. In general,
we recommend the user to start with 2 cores and then adjust accordingly to the
available memory. The code snippet below gives an example of running the grid
search routine over 3 cores.

library(foreach) library(doParallel) ## Fixed control hyperparams mle.control
<\- list(batch.size = 128, epochs = 800, use.GPU = TRUE, early.stopping.epochs
= 50) tune.grid <\- expand.grid(lr = exp(-3:-5), n.knots = c(10,15,20),
n.hidden = c(10,15,20)) set.seed(919) ## Compute CV folds folds <\-
createFolds.SPQR(Y, nfold = 10) cl <\- makeCluster(3) registerDoParallel(cl)
## 10-fold CV for model selection out <\- foreach(i=1:nrow(tune.grid),
.packages = c("SPQR","torch"), .combine = "rbind") %dopar% { lr <\-
tune.grid$lr[i] mle.control$lr <\- lr n.knots <\- tune.grid$n.knots[i]
n.hidden <\- rep(tune.grid$n.hidden[i],2) mle.cv <\- cv.SPQR(X = X, Y = Y,
folds = folds, n.knots = n.knots, n.hidden = n.hidden, method = "MLE", control
= mle.control, normalize = TRUE, verbose = FALSE, seed = 919)
c(lr,n.knots,n.hidden,mle.cv$cve) } stopCluster(cl) ## Print the best model
configuration cve <\- out[,5] best_param <\- out[which.min(cve),]
names(best_param) <\- c("lr","n.knots","n.hidden","n.hidden","cve") best_param
#> lr n.knots n.hidden n.hidden cve #> 6.737947e-03 20 15 15 -161.6610We refit
the model with the best configuration, using 10% of data as validation set for
early stopping criterion.

mle.control$valid.pct <\- 0.1 mle.control$lr <\- exp(-5) mle.fit <\- SPQR(X =
X, Y = Y, n.knots = 20, n.hidden = c(15,15), method = "MLE", control =
mle.control, normalize = TRUE, verbose = TRUE, seed = 919) ## Print model
summary print(mle.fit, showModel=TRUE) #> SPQR fitted using MLE approach #> #>
Learning rate: 0.006737947 #> Batch size: 128 #> #> Model specification: #>
Layers #> Input Output Activation #> 10 15 tanh #> 15 15 tanh #> 15 20 softmax
#> #> Loss: #> train = -219.1049, validation = -169.3304 #> #> Elapsed time:
6.10 minutesTo determine which covariates are important in predicting
quantiles of dem, we compare the quantile VI of all covariates at
$\tau\in\\{0.1,0.5,0.9\\}$. The result is shown in Figure 12.

plotQVI(mle.fit, var.names = colnames(X), tau=c(0.1,0.5,0.9))

Figure 12: Variable importance (VI) of Australia electricity demand data.
Estimated quantile VI of all covariates at $\tau\in\\{0.1,0.5,0.9\\}$.

Across all three quantiles, doy, temp and tod have the highest VI, suggesting
strong seasonal and hourly effects. These effects are also more prominent on
upper tail of the average demand distribution than on central region and lower
tail. While one might expect that there is also a daily/weekend effect on
electric consumption, the estimated effect of dow seems to be almost
negligible, suggesting no particular difference in terms of electric
consumption across the week. One reason could be that only data between 17:30
and 21:30 of the day are used, in which case the probability of people being
home is similar between weekdays and weekends.

To further investigate the effects of doy, temp and tod on the average demand
distribution, we compute their quantile ALE main effects at
$\tau\in\\{0.1,0.5,0.9\\}$. The results are plotted in Figure 13.

for (j in 1:3) plotQALE(mle.fit, var.index=j, tau=c(0.1,0.5,0.9))

Figure 13: Estimated quantile ALE main effects of day of the year (doy), time
of the day (tod) and temperature (temp) at $\tau\in\\{0.1,0.5,0.9\\}$.

The estimated effect of doy has a unimodal shape at $\tau=0.1$ and $\tau=0.5$
and peaks at austral winter. Suggesting a high electric consumption when the
weather is at its coldest. For $\tau=0.9$, the effect also has a minor mode at
austral summer, suggesting that the average demand distribution is highly
right-skewed during austral summer. The estimated effect of tod displays a
quadratic upwarding trend before 20:00 and a linear downwarding trend
afterwards. This suggests that use of electricity becomes more active as
people start to get home and becomes less active as the night closes in. The
downwarding trend is also more prominent at upper tail of the average demand
distribution than at central region or lower tail. This suggests that as the
day ends and people start to go to bed, the probability of very high electric
consumption significantly decreases. The estimated effect of temp seems to
complement that of doy. It has a "check" shape and shows a significant
upwarding trend when temperature is above 20 degree. This suggests a
significantly higher electric consumption when cooling becomes necessary
during summer. We also notice that the turning point for temp effect at upper
quantiles is closer to 15 degree than at lower quantiles. This may suggest
that people who contribute to very high electric consumption tend to turn on
cooling at a lower temperature than those who contribute to very low electric
consumption. Similar characteristics of these effects were also observed in
Fasiolo et al. (2021).

## Summary

In this article, we present the R package SPQR for fitting the semi-parametric
conditional density and quantile regression models as proposed in Xu and Reich
(2021). The main advantage of these models is their capability of modeling
complex covariate effects on the response distribution, which is absent in
existing parametric quantile regression models. Furthermore, the estimated
distribution and quantile functions are always valid, ensuring sensible
inference even when the sample size is small. The package also provides a
framework for fully Bayesian inference of the fitted models to allow
uncertainty quantification, as well as model agnostic tools to understand the
effects of different covariates on different parts of the distribution when
model transparency is needed. We hope that this package will be a valuable
addition to the existing family of quantile regression tools and be especially
useful for analyzing complex and possibly high-dimensional heteroscedastic
data where flexible density and quantile regression models are suitable.
Future work will aim at improving the scalability of the Bayesian SPQR fitting
framework by enabling sparsity-induced priors and implementing variational
Bayes estimator. We also plan to include functions for more specialized use
case, such as estimation of counterfactual distribution and quantile treatment
effects in a causal inference context.

## References

  * Apley (2018) D. Apley.  _ALEPlot: Accumulated Local Effects (ALE) Plots and Partial Dependence (PD) Plots_ , 2018.  URL https://CRAN.R-project.org/package=ALEPlot.  R package version 1.1.
  * Apley and Zhu (2020) D. W. Apley and J. Zhu.  Visualizing the effects of predictor variables in black box supervised learning models.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 82:1059–1086, 2020.
  * Benoit and Van den Poel (2017) D. F. Benoit and D. Van den Poel.  bayesQR: A bayesian approach to quantile regression.  _Journal of Statistical Software_ , 76(7):1–32, 2017.  doi: 10.18637/jss.v076.i07.
  * Betancourt (2017) M. Betancourt.  A conceptual introduction to Hamiltonian Monte Carlo.  _arXiv preprint arXiv:1701.02434_ , 2017.
  * Cannon (2018) A. J. Cannon.  Non-crossing nonlinear regression quantiles by monotone composite quantile regression neural network, with application to rainfall extremes.  _Stochastic Environmental Research and Risk Assessment_ , 32(11):3207–3225. doi:10.1007/s00477–018–1573–6, 2018.
  * Eddelbuettel et al. (2022a) D. Eddelbuettel, R. Francois, J. Allaire, K. Ushey, Q. Kou, N. Russell, I. Ucar, D. Bates, and J. Chambers.  _Rcpp: Seamless R and C++ Integration_ , 2022a.  URL https://CRAN.R-project.org/package=Rcpp.  R package version 1.0.8.
  * Eddelbuettel et al. (2022b) D. Eddelbuettel, R. Francois, D. Bates, B. Ni, and C. Sanderson.  _RcppArmadillo: ’Rcpp’ Integration for the ’Armadillo’ Templated Linear Algebra Library_ , 2022b.  URL https://CRAN.R-project.org/package=RcppArmadillo.  R package version 0.10.8.1.0.
  * Falbel and Luraschi (2022) D. Falbel and J. Luraschi.  _torch: Tensors and Neural Networks with ’GPU’ Acceleration_ , 2022\.  URL https://CRAN.R-project.org/package=torch.  R package version 0.7.2.
  * Fasiolo et al. (2021) M. Fasiolo, S. N. Wood, M. Zaffran, R. Nedellec, and Y. Goude.  qgam: Bayesian nonparametric quantile regression modeling in R.  _Journal of Statistical Software_ , 100(9):1–31, 2021.  doi: 10.18637/jss.v100.i09.
  * Fortuin et al. (2021) V. Fortuin, A. Garriga-Alonso, F. Wenzel, G. Rätsch, R. Turner, M. van der Wilk, and L. Aitchison.  Bayesian neural network priors revisited.  _arXiv preprint arXiv:2102.06571_ , 2021.
  * Galarza et al. (2022) C. E. Galarza, L. Benites, M. Bourguignon, and V. H. Lachos.  _lqr: Robust Linear Quantile Regression_ , 2022.  URL https://CRAN.R-project.org/package=lqr.  R package version 4.1.
  * Geman and Geman (1984) S. Geman and D. Geman.  Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , PAMI-6(6):721–741, 1984.
  * Greenwell et al. (2018) B. M. Greenwell, B. C. Boehmke, and A. J. McCarthy.  A simple and effective model-based variable importance measure.  _arXiv preprint arXiv:1805.04755_ , 2018.
  * Guo et al. (2021) J. Guo, J. Gabry, B. Goodrich, and S. Weber.  _rstan: R Interface to Stan_ , 2021.  URL https://CRAN.R-project.org/package=rstan.  R package version 2.21.3.
  * Hoffman and Gelman (2014) M. D. Hoffman and A. Gelman.  The No-U-Turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo.  _Journal of Machine Learning Research_ , 15:1593–1623, 2014\.
  * Kingma and Ba (2014) D. P. Kingma and J. Ba.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * Koenker (2022) R. Koenker.  _quantreg: Quantile Regression_ , 2022.  URL https://CRAN.R-project.org/package=quantreg.  R package version 5.88.
  * Koenker and Bassett Jr (1978) R. Koenker and G. Bassett Jr.  Regression quantiles.  _Econometrica_ , 46:33–50, 1978.
  * Kuhn (2022) M. Kuhn.  _caret: Classification and Regression Training_ , 2022.  URL https://CRAN.R-project.org/package=caret.  R package version 6.0-91.
  * Lipsitz et al. (2016) M. Lipsitz, A. Belloni, V. Chernozhukov, and I. Fernandez-Val.  _quantreg.nonpar: Nonparametric Series Quantile Regression_ , 2016\.  URL https://CRAN.R-project.org/package=quantreg.nonpar.  R package version 1.0.
  * MacKay (1992) D. J. MacKay.  A Practical Bayesian Framework for Backpropagation Networks.  _Neural Computation_ , 4(3):448–472, 1992.
  * Maidman (2017) A. Maidman.  _plaqr: Partially Linear Additive Quantile Regression_ , 2017.  URL https://CRAN.R-project.org/package=plaqr.  R package version 2.0.
  * Matthews et al. (2018) A. G. d. G. Matthews, M. Rowland, J. Hron, R. E. Turner, and Z. Ghahramani.  Gaussian Process Behaviour in Wide Deep Neural Networks.  _arXiv preprint arXiv:1804.11271_ , 2018.
  * Microsoft and Weston (2022) Microsoft and S. Weston.  _foreach: Provides Foreach Looping Construct_ , 2022.  URL https://CRAN.R-project.org/package=foreach.  R package version 1.5.2.
  * Neal (1996) R. M. Neal.  _Bayesian Learning for Neural Networks_.  Springer, 1996.
  * Neal (2011) R. M. Neal.  MCMC using Hamiltonian dynamics.  In S. Brooks, A. Gelman, G. L. Jones, and X. L. Meng, editors, _Handbook of Markov Chain Monte Carlo_ , chapter 5, pages 113–162. Chapman & Hall/CRC, 2011.
  * Paszke et al. (2019) A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Köpf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, and S. Chintala.  Pytorch: An imperative style, high-performance deep learning library, 2019\.  URL https://arxiv.org/abs/1912.01703.
  * Ramsay (1988) J. O. Ramsay.  Monotone regression splines in action.  _Statistical science_ , pages 425–441, 1988.
  * Stan Development Team (2022) Stan Development Team.  _Stan Modeling Language Users Guide and Reference Manual_ , 2022.  URL https://mc-stan.org.  version 2.29.
  * Tokdar and Cunningham (2019) S. Tokdar and E. Cunningham.  _qrjoint: Joint Estimation in Linear Quantile Regression_ , 2019.  URL https://CRAN.R-project.org/package=qrjoint.  R package version 2.0-3.
  * Vehtari et al. (2017) A. Vehtari, A. Gelman, and J. Gabry.  Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC.  _Statistics and Computing_ , 27(5):1413–1432, 2017.
  * Vehtari et al. (2020) A. Vehtari, J. Gabry, M. Magnusson, Y. Yao, P.-C. Bürkner, T. Paananen, and A. Gelman.  _loo: Efficient Leave-One-Out Cross-Validation and WAIC for Bayesian Models_ , 2020.  URL https://CRAN.R-project.org/package=loo.  R package version 2.4.1.
  * Vladimirova et al. (2019) M. Vladimirova, J. Verbeek, P. Mesejo, and J. Arbel.  Understanding Priors in Bayesian Neural Networks at the Unit Level.  In _International Conference on Machine Learning_ , pages 6458–6467, 2019.
  * Wang and Yan (2021) W. Wang and J. Yan.  _splines2: Regression Spline Functions and Classes_ , 2021.  URL https://CRAN.R-project.org/package=splines2.  R package version 0.4.5.
  * Wickham (2016) H. Wickham.  _ggplot2: Elegant Graphics for Data Analysis_.  Springer-Verlag New York, 2016.  ISBN 978-3-319-24277-4.  URL https://ggplot2.tidyverse.org.
  * Xu and Majumder (2022) S. Xu and R. Majumder.  _SPQR: Semi-Parametric Quantile Regression_ , 2022.  URL https://CRAN.R-project.org/package=SPQR.  R package version 0.1.0.
  * Xu and Reich (2021) S. G. Xu and B. J. Reich.  Bayesian nonparametric quantile process regression and estimation of marginal quantile effects.  _Biometrics_ , 2021.
  * Yu (2021) Y. Yu.  _siqr: An R Package for Single-Index Quantile Regression_ , 2021.  URL https://CRAN.R-project.org/package=siqr.  R package version 0.8.1.

## Acknowledgements

This work was supported by grants from the National Science Foundation
(DMS2152887, CBET2151651), the National Institutes of Health (R01ES031651-01),
and the Southeast National Synthesis Wildfire and the United States Geological
Survey’s National Climate Adaptation Science Center (G21AC10045).

_Steven G. Xu
North Carolina State University
Raleigh, NC
USA
<EMAIL_ADDRESS>
_Reetam Majumder
North Carolina State University
Raleigh, NC
USA
<EMAIL_ADDRESS>
_Brian J. Reich
North Carolina State University
Raleigh, NC
USA
<EMAIL_ADDRESS>