text:
Abstract

Quantiles play a fundamental role in statistics. The quantile function defines the distribution of a random variable and, thus, provides a way to describe the data that is specular but equivalent to that given by the corresponding cumulative distribution function. There are many advantages in working with quantiles, starting from their properties. The renewed interest in their usage seen in the last years is due to the theoretical, methodological, and software contributions that have broadened their applicability. This paper presents the R package Qtools, a collection of utilities for unconditional and conditional quantiles.

Introduction

Quantiles have a long history in applied statistics, especially the median. The analysis of astronomical data by Galileo Galilei in 1632 (Hald 2003, 149) and geodic measurements by Roger Boscovich in 1757 (Koenker 2005, 2) are presumably the earliest examples of application of the least absolute deviation ($$L_1$$) estimator in its, respectively, unconditional and conditional forms. The theoretical studies on quantiles of continuous random variables started to appear in the statistical literature of the 20th century. In the case of discrete data, studies have somewhat lagged behind most probably because of the analytical drawbacks surrounding the discontinuities that characterise discrete quantile functions. Some forms of approximation to continuity have been recently proposed to study the large sample behavior of quantile estimators. For example, Ma, Genton, and Parzen (2011) have demonstrated the asymptotic normality of unconditional sample quantiles based on the definition of the mid-distribution function (Parzen 2004). Machado and Santos Silva (2005) proposed inferential approaches to the estimation of conditional quantiles for counts based on data jittering.

Functions implementing quantile methods can be found in common statistical software. A considerable number of R packages that provide such functions are available on the Comprehensive R Archive Network (CRAN). The base package stats contains basic functions to estimate sample quantiles or compute quantiles of common parametric distributions. The quantreg package (Koenker 2013) is arguably a benchmark for distribution-free estimation of linear quantile regression models, as well as the base for other packages which make use of linear programming (LP) algorithms (Koenker and D’Orey 1987; Koenker and Park 1996). Other contributions to the modelling of conditional quantile functions include packages for Bayesian regression, e.g. bayesQR (Benoit et al. 2014) and BSquare (Smith and Reich 2013), and the lqmm package (Geraci and Bottai 2014; Geraci 2014) for random-effects regression.

The focus of this paper is on the R package Qtools, a collection of models and tools for quantile inference. These include commands for

• quantile-based analysis of the location, scale and shape of a distribution;
• transformation-based quantile regression;
• goodness of fit and restricted quantile regression;
• quantile regression for discrete data;
• quantile-based multiple imputation. The emphasis will be put on the first two topics listed above as they represent the main contribution of the package, while a short description of the other topics is given for completeness.

Unconditional quantiles

Definition and estimation of quantiles

Let $$Y$$ be a random variable with cumulative distribution function (CDF) $$F_{Y}$$ and support $$S_{Y}$$. The CDF calculated at $$y \in S_{Y}$$ returns the probability $$F_{Y}(y) \equiv p = \Pr\left(Y \leq y\right)$$. The quantile function (QF) is defined as $$Q(p) = \inf_{y}\{F_{Y}(y) \geq p\}$$, for $$0 < p < 1$$. (Some authors consider $$0 \leq p \leq 1$$. For practical purposes, it is simpler to exclude the endpoints 0 and 1.) When $$F_{Y}$$ is continuous and strictly monotone (hence, $$f_{Y}(y) \equiv F_{Y}'(y) > 0$$ for all $$y \in S_{Y}$$), the quantile function is simply the inverse of $$F_{Y}$$. In other cases, the quantile $$p$$ is defined, by convention, as the smallest value $$y$$ such that $$F_{Y}(y)$$ is at least $$p$$.

Quantiles enjoy a number of properties. An excellent overview is given by Gilchrist (2000). In particular, the Q-tranformation rule (Gilchrist 2000) or equivariance to monotone transformations states that if $$h(\cdot)$$ is a non-decreasing function on $$\mathbb{R}$$, then $$Q_{h(Y)}(p) = h\left\{Q_{Y}(p)\right\}$$. Hence $$Q_{Y}(p) = h^{-1}\left\{Q_{h(Y)}(p)\right\}$$. Note that this property does not generally hold for the expected value.

Sample quantiles for a random variable $$Y$$ can be calculated in a number of ways, depending on how they are defined (Hyndman and Fan 1996). For example, the function quantile in the base package stats provides nine different sample quantile estimators, which are based on the sample order statistics or the inverse of the empirical CDF. These estimators are distribution-free as they do not depend on any parametric assumption about $$F$$ (or $$Q$$).

Let $$Y_{1}, Y_{2}, \ldots, Y_{n}$$ be a sample of $$n$$ independent and identically distributed (iid) observations from the population $$F_{Y}$$. Let $$\xi_p$$ denote the $$p$$th population quantile and $$\hat{\xi}_p$$ the corresponding sample quantile. (The subscripts will be dropped occasionally to ease notation, e.g. $$F$$ will be used in place of $$F_{Y}$$ or $$\xi$$ in place of $$\xi_{p}$$.) In the continuous case, it is well known that $$\sqrt{n}\left(\hat{\xi}_{p} - \xi_{p}\right)$$ is approximately normal with mean zero and variance $$$\label{eq:1} \tag{1} \omega^2 = \frac{p(1-p)}{\{f_{Y}(\xi_p)\}^2}.$$$ A more general result is obtained when the $$Y_{i}$$’s, $$i = 1, \ldots, n$$, are independent but not identically distributed (nid). The density evaluated at the $$p$$th quantile, $$f(\xi_p)$$, is called density-quantile function by Parzen (1979). Its reciprocal, $$s(p) \equiv 1/f(\xi_p)$$, is called sparsity function (Tukey 1965) or quantile-density function (Parzen 1979).

As mentioned previously, the discontinuities of $$F_{Y}$$ when $$Y$$ is discrete represent a mathematical inconvenience. Ma, Genton, and Parzen (2011) derived the asymptotic distribution of the sample mid-quantiles, that is, the sample quantiles based on the mid-distribution function (mid-CDF). The latter is defined as $$F^{mid}_{Y}(y) = F_{Y}(y) - 0.5p_{Y}(y)$$, where $$p_{Y}(y)$$ denotes the probability mass function (Parzen 2004). In particular, they showed that, as $$n$$ becomes large, $$\sqrt{n}\left(\hat{\xi}^{mid}_{p} - \xi_{p}\right)$$ is approximately normal with mean $$0$$. Under iid assumptions, the expression for the sampling variance is similar to that in (1); see Ma, Genton, and Parzen (2011) for details.

The package Qtools provides the functions midecdf and midquantile, which return objects of class “midecdf” or “midquantile”, respectively, containing: the values or the probabilities at which mid-cumulative probabilities or mid-quantiles are calculated ($$x$$), the mid-cumulative probabilities or the mid-quantiles ($$y$$), and the functions that linearly interpolate those coordinates ($$fn$$). An example is shown below using data simulated from a Poisson distribution.

> library("Qtools")
> set.seed(467)
> y <- rpois(1000, 4)
> pmid <- midecdf(y)
> xmid <- midquantile(y, probs = pmid$y) > pmid Empirical mid-ECDF Call: midecdf(x = y) > xmid Empirical mid-ECDF Call: midquantile(x = y, probs = pmid$y)

A confidence interval for sample mid-quantiles can be obtained using confint.midquantile. This function is applied to the output of midquantile and returns an object of class “data.frame” containing sample mid-quantiles, lower and upper bounds of the confidence intervals of a given level ($$95\%$$ by default), along with standard errors as an attribute named ‘stderr’. This is shown below using the sample $$y$$ generated in the previous example.

> xmid <- midquantile(y, probs = 1:3/4)
> x <- confint(xmid, level = 0.95)
> x

midquantile    lower    upper
25%    2.540000 2.416462 2.663538
50%    3.822816 3.693724 3.951907
75%    5.254902 5.072858 5.436946

> attr(x, "stderr")
[1] 0.06295447 0.06578432 0.09276875
> par(mfrow = c(1,2))
> plot(pmid, xlab = "y", ylab = "CDF", jumps = TRUE)
> points(pmid$x, pmid$y, pch = 15)
> plot(xmid, xlab = "p", ylab = "Quantile", jumps = TRUE)
> points(xmid$x, xmid$y, pch = 15)

Finally, a plot method is available for both midecdf and midquantile objects. An illustration is given in Figure 1. The mid-distribution and mid-quantile functions are discrete and their values are marked by filled squares. The piecewise linear functions connecting the filled squares represent continuous versions of the CDF and QF which interpolate between the steps of, respectively, the ordinary CDF and quantile functions. Note that the argument jumps is a logical value indicating whether values at jumps should be marked.

LSS - Location, scale and shape of a distribution

Since the cumulative distribution and quantile functions are two sides of the same coin, the location, scale, and shape (LSS) of a distribution can be examined using one or the other. Well-known quantile-based measures of location and scale are the median and inter-quartile range (IQR), respectively. Similarly, there are also a number of quantile-based measures for skewness and kurtosis (Groeneveld and Meeden 1984; Groeneveld 1998; Jones, Rosco, and Pewsey 2011).

Define the ‘central’ portion of the distribution as that delimited by the quantiles $$Q(p)$$ and $$Q(1-p)$$, $$0 < p < 0.5$$, and define the ‘tail’ portion as that lying outside these quantiles. Let $$IPR(p) = Q(1-p) - Q(p)$$ denote the inter-quantile range at level $$p$$. Building on the results by Horn (1983) and Ruppert (1987), Staudte (2014) considered the following identity: $$$\label{eq:2} \tag{2} \underbrace{\frac{IPR(p)}{IPR(r)}}_\text{kurtosis} = \underbrace{\frac{IPR(p)}{IPR(q)}}_\text{tail-weight} \, \cdot \, \underbrace{\frac{IPR(q)}{IPR(r)}}_\text{peakedness},$$$ where $$0 < p < q < r < 0.5$$. These quantile-based measures of shape are sign, location and scale invariant. As compared to moment-based indices, they are also more robust to outliers and easier to interpret (Groeneveld 1998; Jones, Rosco, and Pewsey 2011).

It is easy to verify that a quantile function can be written as $$$\label{eq:3} \tag{3} Q(p) = \underbrace{Q(0.5)}_\text{median}\,\, + \,\, \frac{1}{2}\underbrace{IPR(0.25)}_\text{IQR} \, \cdot \, \underbrace{\frac{IPR(p)}{IPR(0.25)}}_\text{shape index} \, \cdot \, \bigg(\underbrace{\frac{Q(p) + Q(1-p) - 2Q(0.5)}{IPR(p)}}_\text{skewness index} - 1\bigg).$$$ This identity establishes a relationship between the location (median), scale (IQR) and shape of a distribution. (This identity appears in Gilchrist (2000, 74) with an error of sign. See also Benjamini and Krieger (1996, eq.1).) The quantity $$IPR(p)/IPR(0.25)$$ in (3) is loosely defined as the ‘shape index’ (Gilchrist 2000, 72), although it can be seen as the tail-weight measure given in (2) when $$p < 0.25$$. For symmetric distributions, the contribution of the skewness index vanishes. Note that the skewness index not only is location and scale invariant, but is also bounded between $$-1$$ and $$1$$ (as opposed to the Pearson’s third standardised moment which can be infinite or even undefined). When this index is near the bounds $$-1$$ or $$1$$, then $$Q(1-p) \approx Q(0.5)$$ or $$Q(p) \approx Q(0.5)$$, respectively.

The function qlss provides a quantile-based LSS summary with the indices defined in (3) of either a theoretical or an empirical distribution. It returns an object of class “qlss”, which is a list containing measures of location (median), scale (IQR and IPR), and shape (skewness and shape indices) for each of the probabilities specified in the argument probs (by default, probs = 0.1). The quantile-based LSS summary of the normal distribution is given in the example below for $$p =0.1$$. The argument fun can take any quantile function whose probability argument is named (this is the case for many standard quantile functions in R, e.g. qt, qchisq, qf, etc.).

> qlss(fun = "qnorm", probs = 0.1)

call:
qlss.default(fun = "qnorm", probs = 0.1)

Unconditional Quantile-Based Location, Scale, and Shape

** Location **
Median
[1] 0
** Scale **
Inter-quartile range (IQR)
[1] 1.34898
Inter-quantile range (IPR)
0.1
2.563103
** Shape **
Skewness index
0.1
0
Shape index
0.1
1.900031

An empirical example is now illustrated using the faithful data set, which contains $$272$$ observations on waiting time (minutes) between eruptions and the duration (minutes) of the eruption for the Old Faithful geyser in Yellowstone National Park, Wyoming, USA. Summary statistics are given in Table 1.

Table 1: Minimum, maximum and three quartiles (Q1, Q2, Q3) for waiting time and duration in the Old Faithful Geyser data set.
 Minimum Q1 Q2 Q3 Maximum Wating time 43.0 58.0 76.0 82.0 96.0 Duration 1.6 2.2 4.0 4.5 5.1

Suppose the interest is in describing the distribution of waiting times. The density is plotted in Figure 2, along with the mid-quantile function. The distribution is bimodal with peaks at around $$54$$ and $$80$$ minutes. Note that the arguments of the base function quantile, including the argument type, can be passed on to qlss.

> x <- as.numeric(faithful$eruptions >= 3) > fit <- rq(formula = y ~ x, tau = c(0.1, 0.25, 0.5, 0.75, 0.9)) > fit Call: rq(formula = y ~ x, tau = c(0.1, 0.25, 0.5, 0.75, 0.9)) Coefficients: tau= 0.10 tau= 0.25 tau= 0.50 tau= 0.75 tau= 0.90 (Intercept) 47 50 54 59 63 x 26 26 26 25 25 Degrees of freedom: 272 total; 270 residual From the output above, it is quite evident that the distribution of waiting times is shifted by an approximately constant amount at all considered values of $$p$$. The location-shift hypothesis can be tested by using the Khmaladze test. The null hypothesis is that two distributions, say $$F_{0}$$ and $$F_{1}$$, differ by a pure location shift (Koenker and Xiao 2002), that is $H_{0}: \, F^{-1}_{1}(p) = F^{-1}_{0}(p) + \delta_{0},$ where $$\delta_{0}$$ is the quantile treatment effect, constant over $$p$$. The location–scale-shift specification of the test considers $H_{0}: \, F^{-1}_{1}(p) = \delta_{1}F^{-1}_{0}(p) + \delta_{0}.$ The alternative hypothesis is that the model is more complex than the one specified in the null hypothesis. The Khmaladze test is implemented in quantreg (see ?quantreg::KhmaladzeTest for further details). The critical values of the test and corresponding significance levels (Koenker 2005) are not readily available in the same package. These have been hardcoded in the Qtools function KhmaladzeFormat which can be applied to “KhmaladzeTest” objects. For the Old Faithful Geyser data, the result of the test is not statistically significant at the $$10\%$$ level. > kt <- KhmaladzeTest(formula = y ~ x, taus = seq(.05, .95, by = .01), > KhmaladzeFormat(kt, 0.05) Khmaladze test for the location-shift hypothesis Joint test is not significant at 10% level Test(s) for individual slopes: not significant at 10% level Goodness of fit Distribution-free quantile regression does not require introducing an assumption on the functional form of the error distribution (Koenker and Bassett 1978), but only weaker quantile restrictions (Powell 1994). Comparatively, the linear specification of the conditional quantile function in Equation 4 is a much stronger assumption and thus plays an important role for inferential purposes. The problem of assessing the goodness of fit (GOF) is rather neglected in applications of QR. Although some approaches to GOF have been proposed (Zheng 1998; Koenker and Machado 1999; X. M. He and Zhu 2003; Khmaladze and Koul 2004), there is currently a shortage of software code available to users. The function GOFTest implements a test based on the cusum process of the gradient vector (X. M. He and Zhu 2003). Briefly, the test statistic is given by the largest eigenvalue of $n^{-1}\sum_{i}^{n} \mathbf{R}_{n}(\mathbf{x}_{i})\mathbf{R}^{\top}_{n}(\mathbf{x}_{i})$ where $$\mathbf{R}_{n}(\mathbf{t}) = n^{-1/2} \sum_{j=1}^{n} \psi_{p}(r_{j})\mathbf{x}_{j} I(\mathbf{x}_{j} \leq \mathbf{t})$$ is the residual cusum (RC) process and $$\psi_{p}(r_{j})$$ is the derivative of the loss function $$\kappa_{p}$$ calculated for residual $$r_{j} = y_{j} - \mathbf{x}_{j}^{\top}\boldsymbol \beta(p)$$. The sampling distribution of this test statistic is non-normal (X. M. He and Zhu 2003) and a resampling approach is used to obtain the $$p$$-value under the null hypothesis. An example is provided further below using the New York Air Quality data set, which contains $$111$$ complete observations on daily mean ozone (parts per billion – ppb) and solar radiation (Langleys – Ly). For simplicity, wind speed and maximum daily temperature, also included in the data set, are not analysed here. Suppose that the model of interest is $$$\label{eq:5} \tag{5} Q_{\text{ozone}}(p) = \beta_{0}(p) + \beta_{1}(p) \cdot \text{Solar.R}.$$$ Three conditional quantiles ($$p \in \{0.1,0.5,0.9\}$$) are estimated and plotted using the following code: > dd <- airquality[complete.cases(airquality), ] > dd <- dd[order(dd$Solar.R), ]
> fit.rq <- rq(Ozone ~ Solar.R, tau = c(.1, .5, .9), data = dd)
> x <- seq(min(dd$Solar.R), max(dd$Solar.R), length = 200)
> yhat <- predict(fit.rq, newdata = data.frame(Solar.R = x))
> plot(Ozone ~ Solar.R, data = dd)
> apply(yhat, 2, function(y,x) lines(x,y), x = x)
## Warning: package 'quantreg' was built under R version 4.1.3

As a function of solar radiation, the median of the ozone daily averages increases by $$0.09$$ ppb for each Ly increase in solar radiation (Figure 4). The 90th centile of conditional ozone shows a steeper slope at $$0.39$$ ppb/Ly, about nine times larger than the slope of the conditional $$10$$th centile at $$0.04$$ ppb/Ly.

The RC test applied to the the object fit.rq provides evidence of lack of fit for all quantiles considered, particularly for $$p = 0.1$$ and $$p = 0.5$$. Therefore the straight-line model in Equation 5 for these three conditional quantiles does not seem to be appropriate. The New York Air Quality data set will be analysed again in the next section, where a transformation-based approach to nonlinear modelling is discussed.

> gof.rq <- GOFTest(fit.rq, alpha = 0.05, B = 1000, seed = 987)
> gof.rq

Goodness-of-fit test for quantile regression based on the cusum process
Quantile 0.1: Test statistic = 0.1057; p-value = 0.001
Quantile 0.5: Test statistic = 0.2191; p-value = 0
Quantile 0.9: Test statistic = 0.0457; p-value = 0.018

Transformation models

Complex dynamics may result in nonlinear effects in the relationship between the covariates and the response variable. For instance, in kinesiology, pharmacokinetics, and enzyme kinetics, the study of the dynamics of an agent in a system involves the estimation of nonlinear models; phenomena like human growth, certain disease mechanisms and the effects of harmful environmental substances such as lead and mercury, may show strong nonlinearities over time. In this section, the linear model is abandoned in favor of a more general model of the type $$$\label{eq:6} \tag{6} Q_{Y|X}(p) = g\left\{\mathbf{x}^{\top}\boldsymbol \beta(p)\right\},$$$ for some real-valued function $$g$$. If $$g$$ is nonlinear, the alternative approaches to conditional quantile modelling are

1. nonlinear parametric models—this approach may provide a model with substantive interpretability, possibly parsimonious (in general more parsimonious than polynomials), and valid beyond the observed range of the data. A nonlinear model depends on either prior knowledge of the phenomenon or the introduction of new, strong theory to explain the observed relationship with potential predictive power. Estimation may present challenges;
2. polynomial models and smoothing splines—this approach goes under the label of nonparametric regression, in which the complexity of the model is approximated by a sequence of linear polynomials (a naive global polynomial trend can be considered to be a special case). A nonparametric model need not introducing strong assumptions about the relationship and is essentially data-driven. Estimation is based on linear approximations and, typically, requires the introduction of a penalty term to control the degree of smoothing; and
3. transformation models—a flexible, parsimonious family of parametric transformations is applied to the response seeking to obtain approximate linearity on the transformed scale. The data provide information about the ‘best’ transformation among a family of transformations. Estimation is facilitated by the application of methods for linear models.

The focus of this section is on the third approach. More specifically the functions available in Qtools refer to the methods for transformation-based QR models developed by Powell (1991), Chamberlain (1994), Mu and He (2007), Dehbi, Cortina-Borja, and Geraci (2016) and Geraci and Jones (2015). Examples of approaches to nonlinear QR based on parametric models or splines can be found in Koenker and Park (1996) and Yu and Jones (1998), respectively.

The goal of transformation-based QR is to fit the model $$$\label{eq:7} \tag{7} Q_{h\left(Y;\lambda_{p}\right)}(p) = \mathbf{x}^{\top}\boldsymbol \beta(p).$$$ The assumption is that the transformation $$h$$ is the inverse of $$g$$, $$h(Y; \lambda_{p}) \equiv g^{-1}(Y)$$, so that the $$p$$th quantile function of the transformed response variable is linear. (In practice, it is satisfactory to achieve approximate linearity.) The parameter $$\lambda_{p}$$ is a low-dimensional parameter that gives some flexibility to the shape of the transformation and is estimated from the data. In general, the interest is on predicting $$Q_{Y|X}(p)$$ and estimating the effects of the covariates on $$Q_{Y|X}(p)$$. If $$h$$ is a non-decreasing function on $$\mathbb{R}$$ (as is the case for all transformations considered here), predictions can be easily obtained from (7) by virtue of the equivariance property of quantiles, $$$\label{eq:8} \tag{8} Q_{Y|X}(p) = h^{-1}\left\{\mathbf{x}^{\top}\boldsymbol \beta(p); \lambda_{p}\right\}.$$$ The marginal effect of the $$j$$th covariate $$x_{j}$$ can be obtained by differentiating the quantile function $$Q_{Y|X}(p)$$ with respect to $$x_{j}$$. This can be written as the derivative of the composition $$Q \circ \eta$$, i.e.
$$$\label{eq:9} \tag{9} \frac{\partial Q(p)}{\partial x_{j}} = \frac{\partial Q(p)}{\partial \eta(p)} \cdot \frac{\partial \eta(p)}{\partial x_{j}},$$$ $$\eta(p) = \mathbf{x}^{\top}\boldsymbol \beta(p)$$. Once the estimates $$\hat{\boldsymbol{\beta}}(p)$$ and $$\hat{\lambda}_{p}$$ are obtained, these can be plugged in Equations 8 and 9.

The package Qtools provides several transformation families, namely the Box–Cox (Box and Cox 1964), Aranda-Ordaz (Aranda-Ordaz 1981), and Jones (Jones 2007; Geraci and Jones 2015) transformations. A distinction between these families is made in terms of the support of the response variable to which the transformation is applied and the number of transformation parameters. The Box–Cox model is a one-parameter family of transformations which applies to singly bounded variables, $$y > 0$$. The Aranda-Ordaz symmetric and asymmetric transformations too have one parameter and are used when responses are bounded on the unit interval, $$0 < y < 1$$ (doubly bounded). Geraci and Jones (2015) developed two families of transformations which can be applied to either singly or doubly bounded responses:

• Proposal I transformations – this family has one parameter and it comes in both symmetric and asymmetric forms;
• Proposal II transformations – this family has two parameters, with one parameter modelling the symmetry (or lack thereof) of the transformation.

Originally, Box and Cox (1964) proposed using power transformations to address lack of linearity, homoscedasticity and normality of the residuals in mean regression modelling. () reported that seldom does this transformation fulfil the basic assumptions of linearity, normality and homoscedasticity simultaneously as originally suggested by Box & Cox (1964). The Box-Cox transformation has found more practical utility in the empirical determination of functional relationships in a variety of fields, especially in econometrics’’.

Indeed, the practical utility of power transformations has been long recognised in QR modelling (Powell 1991; Buchinsky 1995; Chamberlain 1994; Mu and He 2007). Model 7 is the Box–Cox QR model if $$$\label{eq:10} \tag{10} h_{BC}\left(Y;\lambda_{p}\right) = \begin{cases} \dfrac{Y^{\lambda_{p}} - 1}{\lambda_{p}} & \text{if \lambda_{p} \neq 0}\\[.5cm] \log Y & \text{if \lambda_{p} = 0}. \end{cases}$$$ Note that when $$\lambda_{p} \neq 0$$, the range of this transformation is not $$\mathbb{R}$$ but the singly bounded interval $$(-1/\lambda_{p},\infty)$$. This implies that the inversion in (8) is defined only for $$\lambda_{p} \mathbf{x}^{\top}\boldsymbol \beta(p) + 1 > 0$$.

The ‘symmetric’ Aranda-Ordaz transformation is given by $$$\label{eq:11} \tag{11} h_{AOs}\left(Y;\lambda_{p}\right) = \begin{cases} \dfrac{2}{\lambda_{p}} \quad \dfrac{Y^{\lambda_{p}} - \left(1-Y\right)^{\lambda_{p}}}{Y^{\lambda_{p}} + \left(1-Y\right)^{\lambda_{p}}}& \text{if \lambda_{p} \neq 0},\\[.5cm] \log\left(\dfrac{Y}{1-Y}\right) & \text{if \lambda_{p} = 0}. \end{cases}$$$ (The symmetry here is that $$h_{AOs}(\theta;\lambda_p) = -h_{AOs}(1-\theta;\lambda_p)=h_{AOs}(\theta;-\lambda_p)$$.) There is a range problem with this transformation too since, for all $$\lambda_p \neq 0$$, the range of $$h_{AOs}$$ is not $$\mathbb{R}$$, but $$(-2/|\lambda_{p}|, 2/|\lambda_{p}|)$$. The ‘asymmetric’ Aranda-Ordaz transformation is given by $$$\label{eq:12} \tag{12} h_{AOa}\left(Y;\lambda_{p}\right) = \begin{cases} \log \left\{\dfrac{\left(1-Y\right)^{-\lambda_{p}} - 1}{\lambda_{p}} \right\}& \text{if \lambda_{p} \neq 0},\\[.5cm] \log \left\{-\log\left(1 - Y\right)\right\} & \text{if \lambda_{p} = 0}. \end{cases}$$$ For $$\lambda_{p} = 0$$, this is equivalent to the complementary log-log. The asymmetric Aranda-Ordaz transformation does have range $$\mathbb{R}$$. Note that $$h_{AOa}(Y;1) = \log (Y/(1-Y))$$, i.e. the transformation is symmetric.

To overcome range problems, which give rise to computational difficulties, Geraci and Jones (2015) proposed to use instead one-parameter transformations with range $$\mathbb{R}$$. Proposal I is written in terms of the variable (say) $$W$$, where $$$\label{eq:13} \tag{13} h_I\left(W;\lambda_{p}\right) = \begin{cases} \dfrac{1}{2\lambda_{p}}\left(W^{\lambda_{p}} - \dfrac{1}{W^{\lambda_{p}}}\right)& \text{if \lambda_{p} \neq 0}\\[.5cm] \log W & \text{if \lambda_{p} = 0}, \end{cases}$$$ which takes on four forms depending on the relationship of $$W$$ to $$Y$$, as described in Table 2. For each of domains $$(0,\infty)$$ and $$(0,1)$$, there are symmetric and asymmetric forms.

Table 2: Choices of $$W$$ and corresponding notation for transformations based on (13).
 Support of $$Y$$ Symmetric Asymmetric $$(0,\infty)$$ $$W= Y$$ $$W= \log(1+Y)$$ $$h_{Is}(Y;\lambda_p)$$ $$h_{Ia}(Y;\lambda_p)$$ $$(0,1)$$ $$W= Y/(1-Y)$$ $$W= -\log(1-Y)$$ $$h_{Is}(Y;\lambda_p)$$ $$h_{Ia}(Y;\lambda_p)$$

Since the transformation in (13) has range $$\mathbb{R}$$ for all $$\lambda_{p}$$, it admits an explicit inverse transformation. In addition, in the case of a single covariate, every estimated quantile that results will be monotone increasing, decreasing or constant, although different estimated quantiles can have different shapes from this collection. Geraci and Jones (2015) also proposed a transformation that unifies the symmetric and asymmetric versions of $$h_{I}$$ into a single two-parameter transformation, namely $$$\label{eq:14} \tag{14} h_{II}(W;\lambda_p) = h_I(W_{\delta_p};\lambda_p),$$$ where $$h_I$$ is given in (13) and $\begin{equation*} W_{\delta_{p}} = h_{BC}(1+W;\delta_p) = \begin{cases} \dfrac{(1+W)^{\delta_{p}} - 1}{\delta_{p}} & \text{if \delta_{p} > 0}\\[.5cm] \log (1+W) & \text{if \delta_{p} = 0}, \end{cases} \end{equation*}$ with $$W=Y$$, if $$Y > 0$$, and $$W=Y/(1-Y)$$, if $$Y \in (0,1)$$. The additional parameter $$\delta_{p}$$ controls the asymmetry: symmetric forms of $$h_{I}$$ correspond to $$\delta_p = 1$$ while asymmetric forms of $$h_{I}$$ to $$\delta_p = 0$$.

All transformation models discussed above can be fitted using a two-stage (TS) estimator (Chamberlain 1994; Buchinsky 1995) whereby $$\boldsymbol \beta(p)$$ is estimated conditionally on a fine grid of values for the transformation parameter(s). Alternatively, point estimation can be approached using the RC process (Mu and He 2007), which is akin to the process that leads to the RC test introduced in the previous section. The RC estimator avoids the troublesome inversion of the Box-Cox and Aranda-Ordaz transformations, but it is computationally more intensive than the TS estimator.

There are several methods for interval estimation, including those based on large-$$n$$ approximations and the ubiquitous bootstrap. Both the TS and RC estimators have an asymptotic normal distribution. The large-sample properties of the TS estimator for monotonic quantile regression models have been studied by Powell (1991) (see also Chamberlain (1994);Machado and Mata (2000)). Under regularity conditions, it can be shown that the TS estimator is unbiased and will converge to a normal distribution with a sandwich-type limiting covariance matrix which is easy to calculate. In contrast, the form of the covariance matrix of the sampling distribution for the RC estimator is rather complicated and its estimation requires resampling (Mu and He 2007). Finally, if the transformation parameter is assumed to be known, then conditional inference is apposite. In this case, the estimation procedures simplify to those for standard quantile regression problems.

In Qtools, model fitting for one-parameter transformation models can be carried out using the function tsrq. The formula argument specifies the model for the linear predictor as in (7), while the argument tsf provides the desired transformation $$h$$ as specified in Equations 10-13: bc for the Box–Cox model, ao for Aranda-Ordaz families, and mcjI for proposal I transformations. Additional arguments in the function tsrq include

• symmetry, a logical flag to specify the symmetric or asymmetric version of ao and mcjI;
• dbounded, a logical flag to specify whether the response variable is doubly bounded (default is strictly positive, i.e. singly bounded);
• lambda, a numerical vector to define the grid of values for estimating $$\lambda_p$$; and conditional, a logical flag indicating whether $$\lambda_{p}$$ is assumed to be known (in which case the argument lambda provides such known value). \end{itemize}

There are other functions to fit transformation models. The function rcrq fits one-parameter transformation models using the RC estimator. The functions tsrq2 and nlrq2 are specific to Geraci and Jones’s (2015) Proposal II transformations. The former employs a two-way grid search while the latter is based on Nelder-Mead optimization as implemented in optim. Simulation studies in Geraci and Jones (2015) suggest that, although computationally slower, a two grid search is numerically more stable than the derivative-free approach.

A summary of the basic differences between all fitting functions is given in Table 3. The table also shows the available methods in summary.rqt to estimate standard errors and confidence intervals for the model’s parameters. Unconditional inference is carried out jointly on $$\boldsymbol \beta(p)$$ and the transformation parameter by means of bootstrap using the package boot (Canty and Ripley 2014; Davison and Hinkley 1997). Large-$$n$$ approximations (Powell 1991; Chamberlain 1994; Machado and Mata 2000) are also available for the one-parameter TS estimator under iid or nid assumptions.

When summary.rqt is executed with the argument conditional = TRUE, confidence interval estimation for $$\boldsymbol \beta_{p}$$ is performed with one of the several methods developed for linear quantile regression estimators (see options rank, iid, nid, ker, and boot in quantreg::summary.rq).

 Function name Transformation parameters Estimation Standard errors or confidence intervals Unconditional Conditional tsrq 1 Two-stage iid, nid, boot All types rcrq 1 Residual cusum process boot All types tsrq2 2 Two-stage boot All types nlrq2 2 Nelder–Mead boot – Table 3: Transformation-based quantile regression in package Qtools. ‘All types’ consists of options rank, iid, nid, ker, and boot as provided by function summary in package quantreg.

In the New York Air Quality data example, a linear model was found unsuitable to describe the relationship between ozone and solar radiation. At closer inspection, Figure 4 reveals that the conditional distribution of ozone may in fact be nonlinearly associated with solar radiation, at least for some of the conditional quantiles. The model $$$\label{eq:15} \tag{15} Q_{h_{Is}\{\text{ozone}\}}(p) = \beta_{0}(p) + \beta_{1}(p) \cdot \text{Solar.R},$$$ where $$h_{Is}$$ denotes the symmetric version of (13) for a singly bounded response variable, is fitted for the quantiles $$p \in \{0.1,0.5,0.9\}$$ using the following code:

> system.time(fit.rqt <- tsrq(Ozone ~ Solar.R, data = dd, tsf = "mcjI",
+    symm = TRUE, dbounded = FALSE, lambda = seq(1, 3, by = 0.005),
+    conditional = FALSE, tau = c(.1, .5, .9)))
user  system elapsed
0.5     0.0     0.5
> fit.rqt

call:
tsrq(formula = Ozone ~ Solar.R, data = dd, tsf = "mcjI", symm = TRUE,
dbounded = FALSE, lambda = seq(1, 3, by = 0.005), conditional = FALSE,
tau = c(0.1, 0.5, 0.9))

Proposal I symmetric transformation (singly bounded response)

Optimal transformation parameter:
tau = 0.1 tau = 0.5 tau = 0.9
2.210     2.475     1.500

Coefficients linear model (transformed scale):
tau = 0.1  tau = 0.5 tau = 0.9
(Intercept) -3.3357578 -48.737341 16.557327
Solar.R      0.4169697   6.092168  1.443407

Degrees of freedom: 111 total; 109 residual

The TS estimator makes a search for $$\lambda_{p}$$ over the grid $$1.000, 1.005, \ldots, 2.995, 3.000$$. The choice of the search interval usually results from a compromise between accuracy and performance: the coarser the grid, the faster the computation but the less accurate the estimate. A reasonable approach would be to start with a coarse, wide-ranging grid (e.g. seq(-5, 5, by = 0.5)), then center the interval about the resulting estimate using a finer grid, and re-fit the model.

The output above reports the estimates $$\boldsymbol{\hat{\beta}}(p)$$ and $$\hat{\lambda}_p$$ for each quantile level specified in tau. Here, the quantities of interest are the predictions on the ozone scale and the marginal effect of solar radiation, which can obtained using the function predict.rqt.

> x <- seq(9, 334, length = 200)
> qhat <- predict(fit.rqt, newdata = data.frame(Solar.R = x),
+    type = "response")
> dqhat <- predict(fit.rqt, newdata = data.frame(Solar.R = x),
+    type = "maref", namevec = "Solar.R")
The linear component of the marginal effect is calculated as derivative of
Ozone ~ beta1 * Solar.R
with respect to Solar.R

The calculations above are based on a sequence of 200 ozone values in the interval $$[9,334]$$ Ly, as provided via the argument newdata (if this argument is missing, the function returns the fitted values). There are three types of predictions available:

1. link—predictions of conditional quantiles on the transformed scale (7), i.e. $\hat{Q}_{h\left(Y;\hat{\lambda}_{p}\right)}(p) = \mathbf{x}^{T}\hat{\boldsymbol{\beta}}(p);$
2. response—predictions of conditional quantiles on the original scale (8), i.e. $\hat{Q}_{Y|X}(p) = h^{-1}\left\{\mathbf{x}^{\top}\hat{\boldsymbol{\beta}}(p); \hat{\lambda}_{p}\right\};$ and
3. maref—predictions of the marginal effect (9).

In the latter case, the argument namevec is used to specify the name of the covariate with respect to which the marginal effect has to be calculated. The function maref.rqt computes derivatives symbolically using the stats function deriv and these are subsequently evaluated numerically. While the nonlinear component of the marginal effect in Equation 9 (i.e. $$\partial Q(p)/\partial \eta(p)$$) is rather straightforward to derive for any of the transformations (10)-(13), the derivative of the linear predictor (i.e. $$\partial \eta(p)/\partial x_{j}$$) requires parsing the formula argument in order to obtain an expression suitable for deriv. The function maref.rqt can handle simple expressions with common functions like log, exp, etc., interaction terms, and `as is’ terms (i.e. I()). However, using functions that are not recognised by deriv will trigger an error.

The predicted quantiles of ozone and the marginal effects of solar radiation are plotted in Figure 5 using the following code:

> par(mfrow = c(1, 2))
> plot(Ozone ~ Solar.R, data = dd, xlab = "Solar radiation (lang)",
+    ylab = "Ozone (ppb)")
> for(i in 1:3) lines(x, qhat[ ,i], lty = c(1, 2, 4)[i], lwd = 2)
> plot(range(x), range(dqhat), type = "n", xlab = "Solar radiation (lang)",
+    ylab = "Marginal effect")
> for(i in 1:3) lines(x, dqhat[ ,i], lty = c(1, 2, 4)[i], lwd = 2)

The effect of solar radiation on different quantiles of ozone levels shows a nonlinear behavior, especially at lower ranges of radiation (below $$50$$ Ly) and on the median ozone. It might be worth testing the goodness-of-fit of the model. In the previous analysis, it was found evidence of lack of fit for the linear specification (5). In contrast, the output reported below indicates that, in general, the goodness of fit of the quantile models based on the transformation model (15) has improved since the test statistics are now smaller at all values of $$p$$. However, such improvement is not yet satisfactory for the median.

> GOFTest(fit.rqt, alpha = 0.05, B = 1000, seed = 416)

Goodness-of-fit test for quantile regression based on the cusum process
Quantile 0.1: Test statistic = 0.0393; p-value = 0.025
Quantile 0.5: Test statistic = 0.1465; p-value = 0.005
Quantile 0.9: Test statistic = 0.0212; p-value = 0.127

The TS and RC estimators generally provide similar estimates and predictions. However, computation based on the cusum process tends to be somewhat slow, as shown further below. This is also true for the RC test provided by GOFTest.

> system.time(fit.rqt <- rcrq(Ozone ~ Solar.R, data = dd, tsf = "mcjI",
+    symm = TRUE, dbounded = FALSE, lambda = seq(1, 3, by = 0.005),
+    tau = c(.1, .5, .9)))
user  system elapsed
36.88    0.03   37.64

An example using doubly bounded transformations is demonstrated using the A-level Chemistry Scores data set. The latter is available from Qtools and it consists of 31022 observations of A-level scores in Chemistry for England and Wales students, 1997. The data set also includes information of prior academic achievement as assessed with General Certificate of Secondary Education (GCSE) average scores. The goal is to evaluate the ability of GCSE to predict A-level scores. The latter are based on national exams in specific subjects (e.g. chemistry) with grades ranging from A to F. For practical purposes, scores are converted numerically as follows: A = 10, B = 8, C = 6, D = 4, E = 2, and F = 0. The response is therefore doubly bounded between 0 ad 10. It should be noted that this variable is discrete, although, for the sake of simplicity, here it is assumed that the underlying process is continuous.

The model considered here is $$$\label{eq:16} \tag{16} Q_{h_{AOa}\{\text{score}\}}(p) = \beta_{0}(p) + \beta_{1}(p) \cdot \text{gcse},$$$ where $$h_{AOa}$$ denotes the asymmetric Aranda-Ordaz transformation in (12). This model is fitted for $$p = 0.9$$:

> data(Chemistry)
> fit.rqt <- tsrq(score ~ gcse, data = Chemistry, tsf = "ao", symm = FALSE,
+    lambda = seq(0,2,by=0.01), tau = 0.9)

The predicted $$90$$th centile of A-level scores conditional on GCSE and the marginal effect of GCSE are plotted in Figure 6. There is clearly a positive, nonlinear association between the two scores. The nonlinearity is partly explained by the floor and ceiling effects which result from the boundedness of the measurement scale. Note, however, that the S-shaped curve is not symmetric about the inflection point. As a consequence, the marginal effect is skewed to the left. Indeed, the estimate $$\hat{\lambda}_{0.9} = 0$$ and the narrow confidence interval give support to a complementary log-log transformation:

> summary(fit.rqt, conditional = FALSE, se = "nid")

call:
summary.rqt(object = fit.rqt, se = "nid", conditional = FALSE)

Aranda-Ordaz asymmetric transformation (doubly bounded response)

Summary for unconditional inference

tau =  0.9

Optimal transformation parameter:
Value   Std. Error  Lower bound  Upper bound
0.000000000  0.001364422 -0.002674218  0.002674218

Coefficients linear model (transformed scale):
Value  Std. Error Lower bound Upper bound
(Intercept) -4.3520060 0.015414540  -4.3822179  -4.3217941
gcse         0.8978072 0.002917142   0.8920898   0.9035247

Degrees of freedom: 31022 total; 31020 residual

Alternatively, one can estimate the parameter $$\delta_{p}$$ using a two-parameter transformation:

> coef(tsrq2(score ~ gcse, data = chemsub, dbounded = TRUE,
+    lambda = seq(0, 2, by = 0.1), delta = seq(0, 2, by = 0.1),
+    tau = 0.9), all = TRUE)

(Intercept)        gcse      lambda       delta
-4.1442274   0.8681246   0.0000000   0.0000000

These results confirm the asymmetric nature of the relationship since $$\hat{\delta}_{0.9} = 0$$. Similar results (not shown) were obtained with nlrq2.

In conclusion, the package Qtools offers several options in terms of transformations and estimation algorithms, the advantages and disadvantages of which are discussed by . In particular, they found that the symmetric Proposal I transformation improves considerably on the Box-Cox method and marginally on the Aranda-Ordaz transformation in terms of mean squared error of the predictions. Also, asymmetric transformations do not seem to improve sufficiently often on symmetric transformations to be especially recommendable. However, the Box-Cox and the symmetric Aranda-Ordaz transformations should not be used when individual out-of-range predictions represent a potential inconvenience as, for example, in multiple imputation (see section further below). Finally, in some situations transformation-based quantile regression may be competitive as compared to methods based on smoothing, as demonstrated by a recent application to anthropometric charts (Boghossian et al. 2016).

Conditional LSS

Quantile-based measures of location, scale, and shape can be assessed conditionally on covariates. A simple approach is to a fit a linear model as in (4) or a transformation-based model as in (7), and then predict $$\hat{Q}_{Y|X}(p)$$ to obtain the conditional LSS measures in Equation 3 for specific values of $$\mathbf{x}$$.

Estimation of conditional LSS can be carried out by using the function qlss.formula. The conditional model is specified in the argument formula, while the probability $$p$$ is given in probs. (As seen in Equation 3, the other probabilities of interest to obtain the decomposition of the conditional quantiles are $$1-p$$, $$0.25$$, $$0.5$$, and $$0.75$$.) The argument type specifies the required type of regression model, more specifically rq for linear models and rqt for transformation-based models. The function qlss.formula will take any additional argument to be passed to quantreg::rq or tsrq (e.g. subset, weights, etc.).

Let’s consider the New York Air Quality data example discussed in the previous section and assume that the transformation model (15) holds for the quantiles $$p \in \{0.05, 0.1, 0.25, 0.5, 0.75,$$ $$0.9, 0.95\}$$. Then the conditional LSS summary of the distribution of ozone conditional on solar radiation for $$p = 0.05$$ and $$p = 0.1$$ is calculated as follows:

> fit.qlss <- qlss(formula = Ozone ~ Solar.R, data = airquality, type =
+    "rqt", tsf = "mcjI", symm = TRUE, dbounded = FALSE,  lambda =
+    seq(1, 3, by = 0.005), probs = c(0.05, 0.1))
> fit.qlss

call:
qlss.formula(formula = Ozone ~ Solar.R, probs = c(0.05, 0.1),
data = airquality, type = "rqt", tsf = "mcjI", symm = TRUE,
dbounded = FALSE, lambda = seq(1, 3, by = 0.005))

Conditional Quantile-Based Location, Scale, and Shape
-- Values are averaged over observations --

** Location **
Median
[1] 30.2258
** Scale **
Inter-quartile range (IQR)
[1] 43.40648
Inter-quantile range (IPR)
0.05      0.1
88.02909 73.93430
**Shape**
Skewness index
0.05       0.1
0.5497365 0.5180108
Shape index
0.05      0.1
1.960315 1.661648

The output, which is of class “qlss”, is a named list with the same LSS measures seen in the case of unconditional quantiles. However, these are now conditional on solar radiation. By default, the predictions are the fitted values, which are averaged over observations for printing purposes. An optional data frame for predictions can be given via the argument newdata in predict.qlss. If interval = TRUE, the latter computes confidence intervals at the specified level using R bootstrap replications (it is, therefore, advisable to set the seed before calling predict.qlss). The conditional LSS measures can be conveniently plotted using the plot.qlss function as shown in the code below. The argument z is required and specifies the covariate used for plotting. Finally, the argument whichp specifies one probability (and one only) among those given in probs that should be used for plotting (e.g. $$p = 0.1$$ in the following example).

> set.seed(567)
> x <- seq(9, 334, length = 200)
> qhat <- predict(fit.qlss, newdata = data.frame(Solar.R = x),
+    interval = TRUE, level = 0.90, R = 500)
> plot(qhat, z = x, whichp = 0.1, interval = TRUE, type = "l",
+    xlab = "Solar radiation (lang)", lwd = 2)

Figure 7 shows that both the median and the IQR of ozone increase nonlinearly with increasing solar radiation. The distribution of ozone is skewed to the right and the degree of asymmetry is highest at low values of solar radiation. This is due to the extreme curvature of the median which takes on values close to the 10th centile (Figure 5). (Recall that the index approaches 1 when $$Q(p) \approx Q(0.5)$$.) However, the sparsity of observations at the lower end of the observed range of solar radiation determines substantial uncertainty as reflected by the wider confidence interval (Figure 7). At $$p = 0.1$$, the conditional shape index is on average equal to $$1.66$$ and it increases monotonically from $$1.32$$ to about $$1.85$$, remaining always below the tail-weight threshold of a normal distribution ($$1.90$$).

Other functions in Qtools

Restricted quantile regression

Besides a loss of precision, high sparsity (low density) might also lead to a violation of the basic property of monotonicity of quantile functions. Quantile crossing occurs when $$\mathbf{x}_{i}^{\top}\hat{\boldsymbol \beta}(p) > \mathbf{x}_{i}^{\top}\hat{\boldsymbol \beta}(p')$$ for some $$\mathbf{x}_{i}$$ and $$p < p'$$. This problem typically occurs in the outlying regions of the design space (Koenker 2005) where also sparsity occurs more frequently. Balanced designs with larger sample sizes would then offer some assurance against quantile crossing, provided, of course, that the QR models are correctly specified. Model’s misspecification, indeed, can still be a cause of crossing of the quantile curves. Restricted regression quantiles (RRQ) (X. He 1997) might offer a practical solution when little can be done in terms of modelling. This approach applies to a subclass of linear models $Y = \mathbf{x}^{\top}\beta + \epsilon$ and linear heteroscedastic models $Y = \mathbf{x}^{\top}\beta + (\mathbf{x}^{\top}\gamma)\,\epsilon,$ where $$\mathbf{x}^{\top}\gamma > 0$$ and $$\epsilon \sim F$$. Basically, it consists in fitting a reduced regression model passing through the origin. The reader is referred to for details. Here, it is worth stressing that when the restriction does not hold, i.e. if the model is more complex than a location–scale-shift model, then RRQ may yield unsatisfactory results . See also for an examination of the asymptotic properties of the restricted QR estimator. In particular, the relative efficiency of RRQ as compared to RQ depends on the error distribution. For some common unimodal distributions, showed that RRQ in iid models is more efficient than RQ. This property is lost when the error is asymmetric. In contrast, the efficiency of RRQ in heteroscedastic models is comparable to that of RQ even for small samples.

The package Qtools provides the functions rrq, rrq.fit and rrq.wfit which are, respectively, the ‘restricted’ analogous of rq, rq.fit, and rq.wfit in quantreg. S3 methods print, coef, predict, fitted, residuals, and summary are available for objects of class “rrq”. In particular, confidence intervals are obtained using the functions boot and boot.ci from package boot. Future versions of the package will develop the function summary.rrq to include asymptotic standard errors (Zhao 2000). An application is shown below using an example discussed by . The data set, available from Qtools, consists of $$118$$ measurements of esterase concentrations and number of bindings counted in binding experiments.

> data("esterase")
> taus <- c(.1, .25, .5, .75, .9)
> fit.rq <- rq(Count ~ Esterase, data = esterase, tau = taus)
> yhat1 <- fitted(fit.rq)
> fit.rrq <- rrq(Count ~ Esterase, data = esterase, tau = taus)
> yhat2 <- fitted(fit.rrq)

The predicted 90th centile curve crosses the 50th and 75th curves at lower esterase concentrations (Figure 8). The crossing is removed in predictions based on RRQs.

As discussed above, the reliability of the results depends on the validity of the restriction carried by RRQ. A quick check can be performed using the location–scale-shift specification of the Khmaladze test.

> kt <- KhmaladzeTest(formula = Count ~ Esterase, data = esterase,
+ taus = seq(.05,.95,by = .01), nullH = "location-scale")

Khmaladze test for the location-shift hypothesis
Joint test is not significant at 10% level
Test(s) for individual slopes:
not significant at 10% level

The quantile crossing problem can be approached also by directly rearranging the fitted values $$\hat{Q}_{Y|X = \mathbf{x}}(p)$$ to obtain monotone (in $$p$$) predictions for each $$\mathbf{x}$$ (Chernozhukov, Fernandez-Val, and Galichon 2010). This method is implemented in the package Rearrangement (Graybill et al. 2016). As compared to RRQ, this approach is more general as it is not confined to, for example, location–scale-shift models (Chernozhukov, Fernandez-Val, and Galichon 2010); however, in contrast to RRQ, it does not yield estimates of parameters (e.g. slopes) of the model underlying the final monotonised curves. Such estimates, available from rrq objects, may be of practical utility when summarising the results.

Conditional quantiles of discrete data

Modeling discrete response variables, such as categorical and count responses, has been traditionally approached with distribution-based methods: a parametric model $$F_{Y|X}(y;\theta)$$ is assumed and then fitted by means of MLE. Binomial, negative binomial, multinomial and Poisson regressions are well-known in many applied sciences. Because of the computational advantages and the asymptotic properties of MLE, these methods have long ruled among competing alternatives.

Modeling conditional functions of discrete data is less common and, on a superficial level, might even appear as an unnecessary complication. However, a deeper look at its rationale will reveal that a distribution-free analysis can provide insightful information in the discrete case as it does in the continuous case. Indeed, methods for conditional quantiles of continuous distributions can be—and have been—adapted to discrete responses.

Jittering for count responses

Let $$Y$$ be a count variable such as, for example, the number of car accidents during a week or the number of visits of a patient to their doctor during a year. As usual, $$X$$ denotes a vector of covariates. Poisson regression, which belongs to the family of generalized linear models (GLMs), is a common choice for this kind of data, partly because of its availability in many statistical packages. Symbolically, $$Y \sim \textrm{Pois}(\theta)$$ where $\theta \equiv \mathrm{E}(Y|X = \mathbf{x}) = h^{-1}\left(\mathbf{x}^{\top}\boldsymbol{\beta}\right)$ and $$h$$ is the logarithmic link function. Note that the variance also is equal to $$\theta$$. Indeed, moments of order higher than $$2$$ governing the shape of the distribution depend on the same parameter. Every component of the conditional LSS in a Poisson model is therefore controlled by $$\theta$$. If needed, more flexibility can be achieved using a distribution-free approach.

Machado and Santos Silva (2005) proposed the model $$$\label{eq:17} \tag{17} Q_{h\left(Z;p\right)}(p) = \mathbf{x}^{\top}\boldsymbol \beta(p),$$$ where $$Z = Y + U$$ is obtained by jittering $$Y$$ with a $$[0,1)$$-uniform noise $$U$$, independent of $$Y$$ and $$X$$. In principle, any monotone transformation $$h$$ can be considered. A natural choice for count data is a log-linear model , i.e. $\begin{equation*} h\left(Z;p\right) = \begin{cases} \log\left(Z - p\right) & \text{for Z > p}\\[.5cm] \log \zeta & \text{for Z \leq p}. \end{cases} \end{equation*}$ where $$0 < \zeta < p$$. It follows that $$Q_{Z|X}(p) = p + \exp\left(\mathbf{x}^{\top}\boldsymbol \beta(p)\right)$$. (Note that the $$p$$th quantile of the conditional distribution of $$Z$$ is bounded below by $$p$$.) Given the continuity between counts induced by jittering, standard inference for linear quantile functions (Koenker and Bassett 1978) can be applied to fit (17). In practice, a sample of $$M$$ jittered responses $$Z$$ is taken to estimate $$\hat{\boldsymbol \beta}_{m}(p)$$, $$m = 1,\ldots,M$$; the noise is then averaged out, $$\hat{\boldsymbol \beta}(p) = \frac{1}{M} \sum_{m} \hat{\boldsymbol \beta}_{m}(p)$$.

Machado and Santos Silva (2005)’s methods, including large-$$n$$ approximations for standard errors, are implemented in the function rq.counts. The formula argument specifies a log-linear model as in (17). Prior to version 1.4, flexible transformations as described above (e.g., Box-Cox) were allowed. As of version 1.4 of the package, this option has been removed since inference under transformations other than logarithmic was found unreliable and this was likely due to coding issues.

In the example below, estimation is carried out using $$M=50$$ jittered samples and $$\zeta = 10^{-5}$$ (see Machado and Santos Silva (2005) for further details on these settings).

R> data(esterase)
R> fit.rq.counts <- rq.counts(formula = Count ~ Esterase, tau = 0.1,
+ data = esterase, M = 50, zeta = 1e-05)

Figure 9 shows a contrast between centile curves as predicted by the Poisson and the QR models in the Esterase data set. The Poisson distribution clearly underestimates the variability in the data. An empirical modeling of the conditional quantiles seems to be preferred in this case. Of course, the assumption of log-linearity of the models would need to be carefully assessed (note that GOFTest can be applied also to rq.counts objects).

Conditional mid-quantiles for discrete responses

Geraci and Farcomeni (2019) developed quantile regression methods for discrete responses by extending Parzen’s definition of marginal mid-quantiles which we discussed previously. These methods are general as they can be applied to a large variety of discrete responses, including binary, ordinal, and count variables.

Geraci and Farcomeni (2019) first defined the conditional mid-distribution function as $$$\label{eq:18} \tag{18} G_{Y|X}(y) \equiv \Pr(Y \leq y|X) - 0.5\cdot \Pr(Y = y|X).$$$ Then, they introduced the conditional mid-quantile function as its inverse, that is $$H_{Y|X}(p) \equiv G^{-1}_{Y|X}(p)$$.

Estimation is carried out in a two-step approach. In the first step, conditional mid-probabilities are obtained nonparametrically and, in the second step, regression coefficients are estimated by solving an implicit equation. When constraining the quantile index to a data-driven admissible range, the second-step estimating equation has a least-squares type, closed-form solution. Such estimator is shown to be strongly consistent and asymptotically normal. Its good performance is shown in a simulation study with data generated from different discrete response models.

In the example below, the function midrq is applied to the Esterase data set to obtain 10th and 90th conditional mid-quantiles (plotted in Figure 10). As opposed to rq.counts, the function midrq supports Box-Cox and Aranda-Ordaz transformations. S3 methods for midrq objects (e.g., summary, coef, predict, residuals, vcov) are available in the Qtools package. The latter also provides the function cmidecdf to fit the conditional mid-CDF. See Geraci and Farcomeni (2019) for details on estimation.

R> data(esterase)
R> fit.midrq <- midrq(Count ~ Esterase, tau = c(0.5, 0.9),
+ data = esterase, type = 1, lambda = 0)

Quantile-based multiple imputation

Regression models play an important role in conditional imputation of missing values. QR can be used as an effective approach for multiple imputation (MI) when location-shift models are inadequate (Muñoz and Rueda 2009; Bottai and Zhen 2013; Geraci 2016).

In Qtools, mice.impute.rq and mice.impute.rrq are auxiliary functions written to be used along with the functions of the R package mice (van Buuren and Groothuis-Oudshoorn 2011). The former is based on the standard QR estimator (rq.fit) while the latter on the restricted counterpart (rrq.fit). Both imputation functions allow for the specification of the transformation-based QR models discussed previously. The equivariance property is useful to achieve linearity of the conditional model and to ensure that imputations lie within some interval when imputed variables are bounded. An example is available from the help file using the nhanes data set. See also Geraci (2016) and Geraci and McLain (2018) for a thorough description of these methods.

Final remarks

Quantiles have long occupied an important place in statistics. The package Qtools builds on recent methodological and computational developments of quantile functions and related methods to promote their application in statistical data modelling.

Acknowledgements

This work was partially supported by an ASPIRE grant from the Office of the Vice President for Research at the University of South Carolina. Please cite the Qtools package as: Geraci, M. (2016). Qtools: A collection of models and other tools for quantile inference. R Journal, 8(2), 117-138.

References

Aranda-Ordaz, F. J. 1981. “On Two Families of Transformations to Additivity for Binary Response Data.” Biometrika 68 (2): 357–63.
Azzalini, A., and A. W. Bowman. 1990. “A Look at Some Data on the Old Faithful Geyser.” Journal of the Royal Statistical Society C 39 (3): 357–65.
Benjamini, Y., and A. M. Krieger. 1996. “Concepts and Measures for Skewness with Data-Analytic Implications.” Canadian Journal of Statistics 24 (1): 131–40.
Benoit, D. F., R. Al-Hamzawi, K. Yu, and D. Van den Poel. 2014. bayesQR: Bayesian Quantile Regression. https://cran.r-project.org/package=bayesQR.
Boghossian, N. S., M. Geraci, E. M. Edwards, K. A. Morrow, and J. D. Horbar. 2016. “Anthropometric Charts for Infants Born Between 22 and 29 Weeks’ Gestation.” Pediatrics.
Bottai, M., and H. Zhen. 2013. “Multiple Imputation Based on Conditional Quantile Estimation.” Epidemiology, Biostatistics, and Public Health 10 (1): e8758.
Box, G. E. P., and D. R. Cox. 1964. “An Analysis of Transformations.” Journal of the Royal Statistical Society B 26 (2): 211–52.
Buchinsky, M. 1995. “Quantile Regression, Box-Cox Transformation Model, and the US Wage Structure, 1963-1987.” Journal of Econometrics 65 (1): 109–54.
Canty, A., and B. D. Ripley. 2014. Boot: Bootstrap r (s-PLUS) Functions. https://cran.r-project.org/package=boot.
Chamberlain, G. 1994. “Quantile Regression, Censoring, and the Structure of Wages.” In Advances in Econometrics: Sixth World Congress, edited by C. Sims. Vol. 1. Cambridge, UK: Cambridge University Press.
Chernozhukov, V., I. Fernandez-Val, and A. Galichon. 2010. “Quantile and Probability Curves Without Crossing.” Econometrica 78 (3): 1093–1125.
Davison, A. C., and D. V. Hinkley. 1997. Bootstrap Methods and Their Applications. Cambridge: Cambridge University Press.
Dehbi, H.-M., M. Cortina-Borja, and M. Geraci. 2016. “Aranda-Ordaz Quantile Regression for Student Performance Assessment.” Journal of Applied Statistics 43 (1): 58–71.
Doksum, K. 1974. “Empirical Probability Plots and Statistical Inference for Nonlinear Models in the Two-Sample Case.” The Annals of Statistics 2 (2): 267–77.
Geraci, M. 2014. “Linear Quantile Mixed Models: The Lqmm Package for Laplace Quantile Regression.” Journal of Statistical Software 57 (13): 1–29.
———. 2016. “Estimation of Regression Quantiles in Complex Surveys with Data Missing at Random: An Application to Birthweight Determinants.” Statistical Methods in Medical Research 25 (4): 1393–1421.
Geraci, M., and M. Bottai. 2014. “Linear Quantile Mixed Models.” Statistics and Computing 24 (3): 461–79.
Geraci, M., and A. Farcomeni. 2019. “Mid-Quantile Regression for Discrete Responses.” arXiv Preprint arXiv:1907.01945 [Stat.ME].
Geraci, M., and M. C. Jones. 2015. “Improved Transformation-Based Quantile Regression.” Canadian Journal of Statistics 43 (1): 118–32.
Geraci, M., and A. McLain. 2018. “Multiple Imputation for Bounded Variables.” Psychometrika 83 (4): 919–40.
Gilchrist, W. 2000. Statistical Modelling with Quantile Functions. Boca Raton, FL: Chapman & Hall/CRC.
Graybill, W., M. Chen, V. Chernozhukov, I. Fernandez-Val, and A. Galichon. 2016. Rearrangement: Monotonize Point and Interval Functional Estimates by Rearrangement. https://cran.r-project.org/package=Rearrangement.
Groeneveld, R. A. 1998. “A Class of Quantile Measures for Kurtosis.” The American Statistician 52 (4): pp. 325–329.
Groeneveld, R. A., and G. Meeden. 1984. “Measuring Skewness and Kurtosis.” Journal of the Royal Statistical Society D 33 (4): pp. 391–399.
Hald, A. 2003. A History of Probability and Statistics and Their Applications Before 1750. New York, NY: John Wiley & Sons.
He, X. 1997. “Quantile Curves Without Crossing.” The American Statistician 51 (2): 186–92.
He, X. M., and L. X. Zhu. 2003. “A Lack-of-Fit Test for Quantile Regression.” Journal of the American Statistical Association 98 (464): 1013–22.
Horn, P. S. 1983. “A Measure for Peakedness.” The American Statistician 37 (1): 55–56.
Hyndman, R. J., and Y. Fan. 1996. “Sample Quantiles in Statistical Packages.” The American Statistician 50 (4): 361–65.
Jones, M. C. 2007. “Connecting Distributions with Power Tails on the Real Line, the Half Line and the Interval.” International Statistical Review 75 (1): 58–69.
Jones, M. C., J. F. Rosco, and A. Pewsey. 2011. “Skewness-Invariant Measures of Kurtosis.” The American Statistician 65 (2): 89–95.
Khmaladze, E. V., and H. L. Koul. 2004. “Martingale Transforms Goodness-of-Fit Tests in Regression Models,” 995–1034.
Koenker, R. 2005. Quantile Regression. New York, NY: Cambridge University Press.
———. 2013. Quantreg: Quantile Regression. https://cran.r-project.org/package=quantreg.
Koenker, R., and G. Bassett. 1978. “Regression Quantiles.” Econometrica 46 (1): 33–50.
Koenker, R., and V. D’Orey. 1987. “Algorithm AS 229: Computing Regression Quantiles.” Journal of the Royal Statistical Society C 36 (3): 383–93.
Koenker, R., and J. A. F. Machado. 1999. “Goodness of Fit and Related Inference Processes for Quantile Regression.” Journal of the American Statistical Association 94 (448): 1296–1310.
Koenker, R., and B. J. Park. 1996. “An Interior Point Algorithm for Nonlinear Quantile Regression.” Journal of Econometrics 71 (1-2): 265–83.
Koenker, R., and Z. J. Xiao. 2002. “Inference on the Quantile Regression Process.” Econometrica 70 (4): 1583–1612.
Lehmann, E. L. 1975. Nonparametrics: Statistical Methods Based on Ranks. San Francisco, CA: Holden-Day.
Ma, Y., M. G. Genton, and E. Parzen. 2011. “Asymptotic Properties of Sample Quantiles of Discrete Distributions.” Annals of the Institute of Statistical Mathematics 63 (2): 227–43.
Machado, J. A. F., and J. Mata. 2000. “Box–Cox Quantile Regression and the Distribution of Firm Sizes.” Journal of Applied Econometrics 15 (3): 253–74.
Machado, J. A. F., and J. M. C. Santos Silva. 2005. “Quantiles for Counts.” Journal of the American Statistical Association 100 (472): 1226–37.
Mu, Y. M., and X. M. He. 2007. “Power Transformation Toward a Linear Regression Quantile.” Journal of the American Statistical Association 102 (477): 269–79.
Muñoz, J. F., and M. Rueda. 2009. “New Imputation Methods for Missing Data Using Quantiles.” Journal of Computational and Applied Mathematics 232 (2): 305–17.
Parzen, E. 1979. “Nonparametric Statistical Data Modeling.” Journal of the American Statistical Association 74 (365): 105–21.
———. 2004. “Quantile Probability and Statistical Data Modeling.” Statistical Science 19 (4): 652–62.
Powell, J. L. 1991. “Estimation of Monotonic Regression Models Under Quantile Restrictions.” In Nonparametric and Semiparametric Methods in Econometrics and Statistics: Proceedings of the Fifth International Symposium on Economic Theory and Econometrics, edited by W. Barnett, J. Powell, and G. Tauchen, 357–84. New York, NY: Cambridge University Press.
———. 1994. “Estimation of Semiparametric Models.” In Handbook of Econometrics, edited by F. Engle Robert and L. McFadden Daniel, Volume 4:2443–2521. Elsevier.
Ruppert, D. 1987. “What Is Kurtosis?: An Influence Function Approach.” The American Statistician 41 (1): 1–5.
Smith, L., and B. Reich. 2013. BSquare: Bayesian Simultaneous Quantile Regression. https://cran.r-project.org/package=BSquare.
Staudte, R. G. 2014. “Inference for Quantile Measures of Kurtosis, Peakedness and Tail-Weight.” arXiv Preprint arXiv:1047.6461v1 [Math.ST]. https://doi.org/10.1080/03610926.2015.1056366.
Tukey, J. W. 1965. “Which Part of the Sample Contains the Information?” Proceedings of the National Academy of Sciences of the United States of America 53 (1): 127–34.
van Buuren, S., and K. Groothuis-Oudshoorn. 2011. “Mice: Multivariate Imputation by Chained Equations in R.” Journal of Statistical Software 45 (3): 1–67.
Yu, K. M., and M. C. Jones. 1998. “Local Linear Quantile Regression.” Journal of the American Statistical Association 93 (441): 228–37.
Zhao, Q. S. 2000. “Restricted Regression Quantiles.” Journal of Multivariate Analysis 72 (1): 78–99.
Zheng, J. X. 1998. “A Consistent Nonparametric Test of Parametric Regression Models Under Conditional Quantile Restrictions.” Econometric Theory 14 (1): 123–38.


text:
DeepAI

Certified Verification of Relational Properties

The use of function contracts to specify the behavior of functions often remains limited to the scope of a single function call. Relational properties link several function calls together within a single specification. They can express more advanced properties of a given function, such as non-interference, continuity, or monotonicity. They can also relate calls to different functions, for instance, to show that an optimized implementation is equivalent to its original counterpart. However, relational properties cannot be expressed and verified directly in the traditional setting of modular deductive verification. Self-composition has been proposed to overcome this limitation, but it requires complex transformations and additional separation hypotheses for real-life languages with pointers. We propose a novel approach that is not based on code transformation and avoids those drawbacks. It directly applies a verification condition generator to produce logical formulas that must be verified to ensure a given relational property. The approach has been fully formalized and proved sound in the Coq proof assistant.

• 3 publications
• 7 publications
• 8 publications
• 8 publications
01/21/2018

Static and Dynamic Verification of Relational Properties on Self-Composed C Code

Function contracts are a well-established way of formally specifying the...
01/21/2018

Self-composition to Prove Relational Properties in Annotated C Program

Deductive verification provides a powerful tool to show functional prope...
11/26/2018

MetAcsl: Specification and Verification of High-Level Properties

Modular deductive verification is a powerful technique capable to show t...
06/04/2021

Constraint-based Relational Verification

In recent years they have been numerous works that aim to automate relat...
10/20/2019

Relational Test Tables: A Practical Specification Language for Evolution and Security

A wide range of interesting program properties are intrinsically relatio...
07/09/2019

Relational Verification via Invariant-Guided Synchronization

Relational properties describe relationships that hold over multiple exe...
06/03/2022

Prophecy Variables for Hyperproperty Verification

Temporal logics for hyperproperties like HyperLTL use trace quantifiers ...

1 Introduction

Modular deductive verification [DBLP:journals/cacm/Hoare69] allows the user to prove that a function respects its formal specification. More precisely, for a given function , any individual call to can be proved to respect the contract of , that is, basically an implication: if the given precondition is true before the call and the call terminates111Termination can be assumed (partial correctness) or proved separately (full correctness) in a well-known way [Floyd1967]; for the purpose of this paper we can assume it., the given postcondition is true after it. However, some kinds of properties are not easily reducible to a single function call. Indeed, it is frequently necessary to express a property that involves several functions or relates the results of several calls to the same function for different arguments. Such properties are known as relational properties [DBLP:conf/popl/Benton04].

Examples of such relational properties include monotonicity (i.e. ), involving 2 calls, or transitivity (), involving 3 calls. In secure information flow [DBLP:journals/mscs/BartheDR11], non-interference is also a relational property. Namely, given a partition of program variables between high-security variables and low-security variables, a program is said to be non-interferent if any two executions starting from states in which the low-security variables have the same initial values will end up in a final state where the low-security variables have the same values. In other words, high-security variables cannot interfere with low-security ones.

Relational properties can also relate calls to different functions. For instance, in the verification of voting rules [BeckertBormerEA2016], relational properties are used for defining specific properties. Notably, applying the voting rule to a sequence of ballots and a permutation of the same sequence of ballots must lead to the same result, i.e. the order in which the ballots are passed to the voting function should not have any impact on the outcome.

Motivation.

Lack of support for relational properties in verification tools was already faced by industrial users (e.g. in [BishopBC13] for C programs). The usual way to deal with this limitation is to use Self-Composition [DBLP:journals/mscs/BartheDR11, DBLP:conf/fm/SchebenS14, blatterKGP17], Product program [DBLP:conf/fm/BartheCK11] or other Self-Composition optimizations [ShemerCAV2021]. Those techniques are based on code transformations that are relatively tedious and error-prone. Moreover, they are hardly applicable in practice to real-life programs with pointers like in C. Namely, self-composition requires that the compared executions operate on completely separated (i.e. disjoint) memory areas, which might be extremely difficult to ensure for complex programs with pointers.

Example 1 (Motivating Example)

Figure 1 shows an example of two simple C programs performing a swap of the values referred to by pointers x1 and x2 (of type int*). Program uses an auxiliary variable x3 (of type int), while performs an in-place swap using arithmetic operations. As usual in that case, to work correctly, each of these programs needs some separation hypotheses: pointers x1 and x2 should be separated (that is, point to disjoint memory locations) and must not point to x1, x2 themselves and, for , to x3.

Consider a relational property, denoted , stating that both programs, executed from two states in which each of *x1 and *x2 has the same value, will end up in two states also having the same values in these locations. To prove this relational property using self-composition, one typically has to generate a new C program (see Fig. 1) composing and . To avoid name conflicts, we rename their variables by adding, resp., suffixes “_1” and “_2”. The relational property is then expressed by a contract of with a precondition and a postcondition . Obviously, both and must include the equalities: *x1_1==*x1_2 and *x2_1==*x2_2, and must also require the aforementioned separation hypotheses necessary for each function. But for programs with pointers and aliasing, this is not sufficient: the user also has to specify additional separation hypotheses222For convenience of the readers, and are defined in detail in Appendix 0.A. between variables coming from the different programs, that is, in our example, that each of x1_1 and x2_1 is separated from each of x1_2 and x2_2. Without such hypotheses, a deductive verification tool cannot show, for example, that a modification of *x1_1 does not impact *x1_2 in the composed program , and is thus unable to deduce the required property. For real-life programs, such separation hypotheses can be hard to specify or generate. ∎

Approach.

This paper proposes an alternative approach that is not based on code transformation or relational rules. It directly uses a verification condition generator (VCGen) to produce logical formulas to be verified (typically, with an automated prover) to ensure a given relational property. It requires no extra code processing (such as sequential composition of programs or variable renaming). Moreover, no additional separation hypotheses—in addition to those that are anyway needed for each function to work—are required. The locations of each program are separated by construction: each program has its own memory state. The language considered in this work was chosen as a minimal language representative of the main issues relevant for relational property verification: it is a standard While language enriched with annotations, procedures and pointers (see programs and in Fig. 1 for examples). Notably, the presence of dereferencing and address-of operations makes it representative of various aliasing problems with (possibly, multiple) pointer dereferences of a real-life language like C. We formalize the proposed approach and prove333The full development is available at https://bit.ly/3FJhE41. its soundness in the proof assistant [Coq]. Our development contains about 3400 lines.

Contributions.

The contributions of this paper include:

• a formalization and proof of soundness of recursive Hoare triple verification with a verification condition generator on a representative language with procedures and aliasing;

• a novel method for verifying relational properties using a verification condition generator, without relying on code transformation (such as self-composition) or making additional separation hypotheses in case of aliasing;

• a formalization and proof of soundness of the proposed method of relational property verification for the considered language.

Outline.

Section 2 introduces an imperative language used in this work. Functional correctness is defined in Section 3, and relational properties in Section 4. Then, we prove the soundness of a VCGen in Section 5, and show how it can be soundly extended to verify relational properties in Section 6. Finally, we present related work in Section 7 and concluding remarks in Section 8.

2 Syntax and Semantics of the Language

2.1 Notation for Locations, States, and Procedure Contracts

We denote by the set of natural numbers, by the set of nonzero natural numbers, and by the set of Boolean values. Let be the set of program locations and the set of program (procedure) names, and let and denote metavariables ranging over those respective sets. We assume that there exists a bijective function , so that . Intuitively, we can see as the address of location .

Let be the set of functions , called memory states, and let denote metavariables ranging over the set. A state maps a location to a value using its address: location has value

We define the update operation of a memory state , also denoted by , as the memory state mapping each address to the same value as , except for , bound to . Formally, is defined by the following rules:

 ∀\memvar∈\Memvar,\locvali∈\Locval,n∈\Natset,\locvalj∈\Locval. i=j⇒σ[i/n](j)=n, (1) ∀\memvar∈\Memvar,\locvali∈\Locval,n∈\Natset,\locvalj∈\Locval. i≠j⇒σ[i/n](j)=\memvar(j). (2)

Let be the set of functions , called procedure environments, mapping program names to commands (defined below), and let denote metavariables ranging over . We write to refer to , the commands (or body) of procedure for a given procedure environment .

Assertions are predicates of arity one, taking as parameter a memory state and returning an equational first-order logic formula. Let metavariables range over the set of assertions. For instance, using -notation, assertion assessing that location is bound to can be defined by This form will be more convenient for relational properties (than e.g. ) as it makes explicit the memory states on which a property is evaluated.

Finally, we define the set of contract environments , and metavariables to range over . More precisely, maps a procedure name to the associated (procedure) contract , composed of a pre- and a postcondition for procedure . As usual, a procedure contract will allow us to specify the behavior of a single procedure call, that is, if we start executing in a memory state satisfying , and the evaluation terminates, the final state satisfies .

2.2 Syntax for Expressions and Commands

Let , and denote respectively the sets of arithmetic expressions, Boolean expressions and commands. We denote by ; and metavariables ranging, respectively, over those sets. Syntax of arithmetic and Boolean expressions is given in Fig. 2. Constants are natural numbers or Boolean values. Expressions use standard arithmetic, comparison and logic binary operators, denoted respectively , , . Since we use natural values, the subtraction is bounded by 0, as in : if , the result of is considered to be 0. Expressions also include locations, possibly with a dereference or address operators.

Figure 2 also presents the syntax of commands in . Sequences, skip and conditions are standard. An assignment can be done to a location directly or after a dereference. Recall that a location contains as a value a natural number, say , that can be seen in turn as the address of a location, namely , so the assignment writes the value of expression to the location , while the address operation computes the address of . An assertion command indicates that an assertion should be valid at the point where the command occurs. The loop command is always annotated with an invariant . As usual, this invariant should hold when we reach the command and be preserved by each loop step. Command is a procedure call. All annotations (assertions, loop invariants and procedure contracts) will be ignored during the program execution and will be relevant only for program verification in Section 5. Procedures do not have explicit parameters and return values (hence we use the term procedure call rather than function call). Instead, as in assembly code [Irvine:2014:ALX:2655333], parameters and return value(s) are shared implicitly between the caller and the callee through memory locations: the caller must put/read the right values at the right locations before/after the call. Finally, to avoid ambiguity, we regroup sequences of commands with .

Example 2

Figure 3 shows an example of a command and a procedure environment where procedure points to a recursive command, called in . With the semantics of Sec. 2.3, from any initial state, the command will return a state in which . Procedure returns a state where if the initial state satisfies . This can be expressed by the contract environment given (in -notation) in Fig. 3. ∎

2.3 Operational Semantics

Evaluation of arithmetic and Boolean expressions in is defined by functions and . Selected evaluation rules for arithmetic expressions are shown in Fig. 4. Operations and have a semantics similar to the C language, i.e. dereferencing and address-of. Semantics of Boolean expressions is standard [DBLP:books/daglib/0070910].

Based on these evaluation functions, we can define the operational semantics of commands in a given procedure environment . Selected evaluation rules444For convenience of the readers, full versions of Fig. 45 are given in Appendix 0.B. are shown in Fig. 5. As said above, both assertions and loop invariants can be seen as program annotations that do not influence the execution of the program itself. Hence, command is equivalent to a skip. Likewise, loop invariant has no influence on the semantics of .

We write to denote that can be derived from the rules of Fig. 5. Our formalization, inspired by [SF], provides a deep embedding of , with an associated parser, in files Aexp.v, Bexp.v and Com.v.

3 Functional Correctness

We define functional correctness in a similar way to the original Hoare triple definition [DBLP:journals/cacm/Hoare69], except that we also need a procedure environment , leading to a quadruple denoted . We will however still refer by the term “Hoare triple” to the corresponding program property, formally defined as follows.

Definition 1 (Hoare triple)

Let be a command, a procedure environment, and and two assertions. We define a Hoare triple as follows:

 \hoaredP\comQ\memcom≜∀\memvar,\memvar′∈\Memvar. \assertP\memvar∧(⊩\eval\com\memvar\memcom\memvar′)⇒\assertQ\memvar′.

Informally, our definition states that, for a given , if a state satisfies and the execution of on terminates in a state , then satisfies .

Next, we introduce notation to denote the fact that, for the given and , every procedure satisfies its contract.

Definition 2 (Contract Validity)

Let be a procedure environment and a contract environment. We define contract validity as follows:

 \hoarectx\memcom\memann≜∀\loccom∈\Loccom. \hoared\pre\loccom\memann\wcall\loccom\post\loccom\memann\memcom).

The notion of contract validity is at the heart of modular verification, since it allows assuming that the contracts of the callees are satisfied during the verification of a Hoare triple. More precisely, to state the validity of procedure contracts without assuming anything about their bodies in our formalization, we will consider an arbitrary choice of implementations of procedures that satisfy the contracts, like in assumption (3) in Lemma 1. This technical lemma, taken from [DBLP:series/txcs/AptBO09, Equation (4.6)], gives an alternative criterion for validity of procedure contracts: if, under the assumption that the contracts in hold, we can prove for each procedure that its body satisfies its contract, then the contracts are valid.

Given a procedure environment and a contract environment such that

 ∀\memcom′∈\Memcom. \hoaredctx\pre\loccom\memann\body\loccom\memcom\post\loccom\memann∀\loccom∈\Loccom,\memcom′\hoarectx\memcom′\memann, (3)

we have

Proof

Any given terminating execution traverses a finite number of procedure calls (over all procedures) that can be replaced by inlining the bodies a sufficient number of times. We first formalize a theory of -inliners (that inline procedure bodies a finite number of times and replace deeper calls by nonterminating loops) and prove their properties. Relying on this elegant theory, the proof of the lemma proceeds by induction on the number of procedure inlinings. ∎

From that, we can establish the main result of this section. Theorem 3.1, taken from [DBLP:series/txcs/AptBO09, Th. 4.2] states that holds if assumption (3) holds and if the validity of contracts of for implies the Hoare triple itself. This theorem is the basis for modular verification of Hoare Triples, as done for instance in Hoare Logic [DBLP:journals/cacm/Hoare69, DBLP:books/daglib/0070910] or verification condition generation.

Theorem 3.1 (Recursion)

Given a procedure environment and a contract environment such that

 ∀\memcom′∈\Memcom. \hoaredctx\pre\loccom\memann\body\loccom\memcom\post\loccom\memann∀\loccom∈\Loccom,\memcom′\hoarectx\memcom′\memann, and \hoaredctxP\comQ\memcom\hoarectx\memcom\memann,

we have

Proof

By Lemma 1.∎

We refer the reader to the development, more precisely the results recursive_proc and recursive_hoare_triple in file Hoare_Triple.v for complete proofs of Lemma 1 and Theorem 3.1 for . To the best of our knowledge, this is the first mechanized proof of these classical results in .

An interesting corollary can be deduced from Theorem 3.1.

[Procedure Recursion] Given a procedure environment and a contract environment such that

 ∀\memcom′∈\Memcom. \hoaredctx\pre\loccom\memann\body\loccom\memcom\post\loccom\memann∀\loccom∈\Loccom,\memcom′\hoarectx\memcom′\memann,

we have

4 Relational Properties

Relational properties can be seen as an extension of Hoare triples. But, instead of linking one program with two properties, the pre- and postconditions, relational properties link programs to two properties, called relational assertions. We define a relational assertion as a predicate taking a sequence of memory states and returning a first-order logic formula. We use metavariables to range over the set of relational assertions, denoted . As a simple example of a relational assertion, we might say that two states bind location to the same value. This would be stated as follows: .

A relational property is a property about programs , stating that if each program starts in a state and ends in a state such that holds, then holds, where and are relational assertions over memory states.

We formally define relational correctness similarly to functional correctness (cf. Def. 1), except that we now use sequences of memory states and commands of equal size. We denote by a sequence of elements where ranges from to . If , is the empty sequence denoted .

Definition 3 (Relational Correctness)

Let be a procedure environment, a sequence of commands (), and and two relational assertions over states. The relational correctness of with respect to and , denoted , is defined as follows:

 \rhoaredP\sequence\comnQ\memcom≜ ∀\sequence\memvar\nat,\sequence\memvar′\nat. \rassertP\sequence\memvar\nat∧(\evalrr\com\memvar\memvar′ψ\nat)⇒\rassertQ\sequence\memvar′\nat.

This notation generalizes the one proposed by Benton [DBLP:conf/popl/Benton04] for relational properties linking two commands: As Benton’s work mostly focused on comparing equivalent programs, using symbol was quite natural.

Example 3 (Relational property)

Figure 6 formalizes the relational property for programs and discussed in Ex. 1. Recall that (written in Fig. 6 in Benton’s notation) states that both programs executed from two states named and having the same values in and will end up in two states and also having the same values in these locations. Notice that the initial state of each program needs separation hypotheses (cf. the second line of the definition of ). Namely, and must point to different locations and must not point to , or, for , to for the property to hold. This relational property is formalized in the developement in file Examples.v. ∎

5 Verification Condition Generation for Hoare Triples

A standard way [Floyd1967] for verifying that a Hoare triple holds is to use a verification condition generator (VCGen). In this section, we formalize a VCGen for Hoare triples and show that it is correct, in the sense that if all verification conditions that it generates are valid, then the Hoare triple is valid according to Def. 1.

5.1 Verification Condition Generator

We have chosen to split the VCGen in three steps, as it is commonly done [Frama-C]:

• function generates the main verification condition, expressing that the postcondition holds in the final state, assuming auxiliary annotations hold;

• function generates auxiliary verification conditions stemming from assertions, loop invariants, and preconditions of called procedures;

• finally, function generates verification conditions for the auxiliary procedures that are called by the main program, to ensure that their bodies respect their contracts.

Definition 4 (Function \Tcn generating the main verification condition)

Given a command , a memory state representing the state before the command, a contract environment , and an assertion , function returns a formula defined by case analysis on as shown in Fig. 7.

Assertion represents the postcondition we want to verify after the command executed from state . For each command, except sequence and branch, a fresh memory state is introduced and related to the current memory state . The new memory state is given as parameter to . For , which does nothing, both states are identical. For assignments, is simply the update of . An assertion introduces a hypothesis over but leaves it unchanged. For a sequence, we simply compose the conditions, that is, we check that the final state of is such that will be verified after executing . For a conditional, we check that if the condition evaluates to , the then branch will ensure the postcondition, and that otherwise the else branch will ensure the postcondition. The rule for calls simply assumes that verifies . Finally, assumes that, after a loop, is a state where the loop condition is and the loop invariant holds. As for an assertion, the callee’s precondition and the loop invariant are just assumed to be true; function , defined below, generates the corresponding proof obligations.

Example 4

For , and , we have:

 \Tcnf\com\mem\memannf≡∀\mem′1.\mem=\mem′1⇒(∀\mem′2.\mem′2=\setmem\mem′112⇒\mem′2(1)=2).\squareforqed
Definition 5 (Function \Tan generating the auxiliary verification conditions)

Given a command , a memory state representing the state before the command, and a contract environment , function returns a formula defined by case analysis on as shown in Fig. 8.

Basically, collects all assertions, preconditions of called procedures, as well as invariant establishment and preservation, and lifts the corresponding formulas to constraints on the initial state through the use of .

Finally, we define the function for generating the conditions for verifying that the body of each procedure defined in respects its contract defined in .

Definition 6 (Function \Tfn generating the procedure verification condition)

Given two environments and , returns the following formula:

 \Tfnf\memann\memcom≜∀\loccom,\mem.\pre\loccom\memann(\mem)⇒\Tanf\body\loccom\memcom\mem\memann∧\Tcnf\body\loccom\memcom\mem\memann\post\loccom\memann.

The VCGen is defined in file Vcg.v of the development. Interested readers will also find a proof (in file Vcg_Opt.v) of a VCGen optimization (not detailed here), which prevents the size of the generated formulas from becoming exponential in the number of conditions in the program [DBLP:conf/popl/FlanaganS01], which is a classical problem for “naive” VCGens.

5.2 Hoare Triple Verification

We can now state the theorems establishing correctness of the VCGen. Their proof can be found in file Correct.v of the development.

First, Lemma 2 shows that, under the assumption of the procedure contracts, a Hoare triple is valid if for all memory states satisfying the precondition, the main verification condition and the auxiliary verification conditions hold.

Lemma 2

Assume the following two properties hold:

 ∀\mem∈Σ,\assertP\mem⇒\Tanf\com\mem\memann, ∀\mem∈Σ,\assertP\mem⇒\Tcnf\com\mem\memannQ.

Then we have

Proof

By structural induction over . ∎

Next, we prove in Lemma 3 that if holds, then for an arbitrary choice of implementations of procedures respecting the procedure contracts, the body of each procedure respects its contract.

Lemma 3

Assume that the formula is satisfied. Then we have

 ∀\memcom′∈\Memcom. \hoaredctx\pre\loccom\memann\body\loccom\memcom\post\loccom\memann∀\loccom∈\Loccom,\memcom′\hoarectx\memcom′\memann.
Proof

By Lemma 2. ∎

Finally, we can establish the main theorem of this section, stating that the VCGen is correct with respect to our definition of Hoare triples.

Theorem 5.1 (Soundness of VCGen)

Assume that we have and

 ∀\mem∈Σ,\assertP\mem⇒\Tanf\com\mem\memann, ∀\mem∈Σ,\assertP\mem⇒\Tcnf\com\mem\memannQ.

Then we have .

Proof

By Theorem 3.1 and Lemmas 2 and 3. ∎

Example 5

Consider again the command , procedure environment , and contract environment of Ex. 2 (presented in Fig. 3). We can apply Theorem 5.1 to prove its functional correctness expressed by the following Hoare triple:

 \memcom:{λ\mem.True}\crec{λ\mem.\mem(2)=\mem(4)×\mem(3)}

(see command com_rec in file Examples.v). ∎

6 Relational Properties Verification

In this section, we propose a verification method for relational properties (defined in Section 4) using the VCGen defined in Section 5 (or, more generally, any VCGen respecting Theorem 5.1). First, we define the notation for the recursive call of function on a sequence of commands and memory states:

Definition 7 (Function \Tr)

Given a sequence of commands and a sequence of memory states , a contract environment and a predicate over states, function is defined by induction on as follows.

• Basis: .

 \Trf[ ][ ]\memann\relaQ≜\rassertQ[ ].
• Inductive: .

 \Trf\sequence\com\nat\sequence\mem\nat\memann\relaQ≜ \Tcnf\com\nat\mem\nat\memannλ\mem′\nat.\Trf\sequence\com\nat−1\sequence\mem\nat−1\memannλ\sequence\mem′\nat−1.\rassertQ\sequence\mem′\nat.

Intuitively, for , gives the weakest relational condition that and must fulfill in order for to hold after executing from and from :

Remark 1

Assume we have a command , a sequence of commands , and a sequence of memory states . From Def. 1, it follows that

 ∀\mem\nat,\mem′\nat. \rassertP\sequence\mem\nat∧(⊩\eval\com\nat\mem\nat\memcom\mem′\nat)⇒ \Trf\sequence\com\nat−1\sequence\mem\nat−1\memannλ\sequence\mem′\nat−1.\rassertQ\sequence\mem′\nat

is equivalent to

 \hoaredλ\memn.\rassertP\sequence\mem\nat\com\natλ\mem′\nat.\Trf\sequence\com\nat−1\sequence\mem\nat−1\memannλ\sequence\mem′\nat−1.\rassertQ\sequence\mem′\nat\memcom.
Example 6 (Relational verification condition)

In order to make things more concrete, we can go back to the relational property between two implementations and of swap defined in Ex. 1 and examine what would be the main verification condition generated by . Let and be defined as in Ex. 3. In this particular case, we have , and is empty (since we do not have any function call), thus Def. 7 becomes:

 \Trf(\cswI,\cswII)(\mem1,\mem2)∅\relaQ=\Tcnf\cswII\mem2∅λ\mem′2.\Tcnf\cswI\mem1∅λ\mem′1.\rassertQ\mem′1,\mem′2.

We thus start by applying over , to obtain, using the rules of Def. 4 for sequence and assignment, the following intermediate formula:

 \Trf(\cswI,\cswII)(\mem1,\mem2)∅\relaQ= \Tcn(\cswII,\mem2,∅, λ\mem′2.∀\mem3,\mem5,\mem7. \mem3=\mem1[3/\mem1(\mem1(1))]⇒ \mem5=\mem3[\mem3(1)/\mem3(\mem3(2))]⇒ \mem7=\mem5[\mem5(2)/\mem5(3)]⇒\rassertQ\mem7,\mem′2.

We can then do the same with to obtain the final formula:

 \Trf(\cswI,\cswII)(\mem1,\mem2)∅\relaQ= ∀\sequence\mem8. \mem4=\mem2[\mem2(1)/\mem2(\mem2(1))+\mem2(\mem2(2))]⇒ \mem6=\mem4[\mem4(2)/\mem4(\mem4(1))−\mem4(\mem4(2))]⇒ \mem8=\mem6[\mem6(1)/\mem6(\mem6(1))−\mem6(\mem6(2))]⇒ \mem3=\mem1[3/\mem1(\mem1(1))]⇒ \mem5=\mem3[\mem3(1)/\mem3(\mem3(2))]⇒ \mem7=\mem5[\mem5(2)/\mem5(3)]⇒\rassertQ\mem7,\mem8.

Here,

with odd (resp., even) indices result from

for (resp., ). ∎

We similarly define a notation for the auxiliary verification conditions for a sequence of commands.

Definition 8 (Function \Tar)

Given a sequence of commands and a sequence of memory states , we define function as follows:

 \Tarf\sequence\com\nat\sequence\mem\nat\memann≜n⋀i=1\Tanf\comi\memi\memann.
Remark 2

For it trivially follows from Def. 8 that:

 \Tarf\sequence\com\nat\sequence\mem\nat\memann≡\Tanf\com\nat\mem\nat\memann∧\Tarf\sequence\com\nat−1\sequence\mem\nat−1\memann.

Using functions and , we can now give the main result of this paper: it states that the verification of relational properties using the VCGen is correct.

Theorem 6.1 (Soundness of relational VCGen)

For any sequence of commands , contract environment , procedure environment , and relational assertions over states and , if the following three properties hold:

 \Tfnf\memann\memcom, (4)
 ∀\sequence\mem\nat,\rassertP\sequence\mem\nat⇒\Tarf\sequence\com\nat\sequence\mem\nat\memann, (5)
 ∀\sequence\mem\nat,\rassertP\sequence\mem\nat⇒\Trf\sequence\com\nat\sequence\mem\nat\memann\relaQ, (6)

then we have

In other words, a relational property is valid if all procedure contracts are valid, and, assuming the relational precondition holds, both the auxiliary verification conditions and the main relational verification condition hold. We give the main steps of the proof below. The corresponding formalization is available in file Rela.v, and the proof of Theorem 6.1 is in file Correct_Rela.v.

Proof

By induction on the length of the sequence of commands .

• Induction basis: . By Def. 3, our goal becomes:

 \rhoaredP\sequence\com0Q\memcom≡\rassertP[ ]⇒\rassertQ[ ].

Indeed, by definition of and Hypothesis (6), holds.

• Induction step: assuming the result for , we prove it for . So, assume we have a sequence of commands , relational assertions and environments respecting (4), (5), (6) (stated for sequences of elements). We have to prove , which, by Def. 3, is equivalent to:

 ∀\sequence\memvar\nat+1,\sequence\memvar′\nat+1. \rassertP\sequence\memvar\nat+1∧

format: plain


text:
### Sidebar

seminars:sml:160426

Statistical Machine Learning Seminar
Hosted by Department of Mathematical Sciences

• Date: Tuesday, April 26, 2016
• Time: 12:00-1:00
• Room: WH-100E
• Speaker: Wolfgang Wefelmeyer (Universität zu Köln)
• Title: Density estimators in regression models with errors in covariates

Abstract

In regression models $Y=r(X)+\varepsilon$ with $X$ and $\varepsilon$ independent, the density of the response $Y$ can be estimated by a convolution of (kernel) estimators for the densities of $r(X)$ and $\varepsilon$. The rate of this convolution estimator depends on the smoothness of the densities of $X$ and $\varepsilon$ and on the smoothness and local flatness of the regression function $r$. When we observe the covariates $X$ with measurement errors, $Z=X+\eta$, we need deconvolution estimators for the densities of $X$ and $\varepsilon$ and for $r$. This is joint work with Anton Schick and Ursula U. Müller.