# New Link Functions for Distribution–Specific Quantile Regression Based on Vector Generalized Linear and Additive Models

Authors: V. F. Miranda-Soberanis; T. W. Yee
Journal: Journal of Probability and Statistics
 (2019)

---

## Abstract

In the usual quantile regression setting, the distribution of the response given the explanatory variables is unspecified. In this work, the distribution is specified and we introduce new link functions to directly model specified quantiles of seven 1–parameter continuous distributions. Using the vector generalized linear and additive model (VGLM/VGAM) framework, we transform certain prespecified quantiles to become linear or additive predictors. Our parametric quantile regression approach adopts VGLMs/VGAMs because they can handle multiple linear predictors and encompass many distributions beyond the exponential family. Coupled with the ability to fit smoothers, the underlying strong assumption of the distribution can be relaxed so as to offer a semiparametric–type analysis. By allowing multiple linear and additive predictors simultaneously, the quantile crossing problem can be avoided by enforcing parallelism constraint matrices. This article gives details of a software implementation called theVGAMextra package forR. Both the data and recently developed software used in this paper are freely downloadable from the internet.


---

## Body

## 1. Introduction

### 1.1. Background

Much of modern regression analysis for estimating conditional quantile functions may be viewed as starting from Koenker and Bassett [1], who offered a systematic strategy for examining how covariates influence the entire response distribution. The fundamental idea is based on the linear specification of the τth quantile function Qy(τ∣x)=βτTx and finding βτ∈Rp that solves the optimization problem(1)minβτ∈Rp⁡∑ρτyi-βτTxi,for independent and identically distributed (i.i.d.) observations from a family of linear quantile regression models, say yi=βτTx+εi,τ,i=1,…,n. Equation (1) can be reformulated as a linear programming problem using the piecewise linear function ρτ(u)=u·τ-I(u<0) for τ∈(0,1). More details can be found in Koenker [2].In the spirit of quantile regression, the conditional distributionY∣x is usually unspecified, although it relies on normal–based asymptotic theory that is used for inference, whilst the assumption of homoskedasticity of the error terms εi,τ is dropped. In this paper we use an alternative approach of conditional–quantile regression based on assuming a prespecified distribution for the response. Parametric quantile regression has some advantages over many nonparametric approaches, including overcoming the quantile crossing problem. Two examples are Noufaily and Jones [3] which is based on the generalized gamma distribution and generalized additive models for location, scale, and shape (GAMLSS; [4]). Further examples are the LMS-BCN method involving the standard normal distribution and a three–parameter Box–Cox transformation [5] and the classical method of quantile regression based on the asymmetric Laplace distribution (ALD).Our approach uses the vector generalized linear and additive model (VGLM/VGAM; [6, 7]) framework. We develop new link functions, G, for the quantile regression model(2)Y∣x~Fx;θ,(3)ητθ=GQyτ∣x,θ,for a vector of quantiles τ=τ1,…,τLT. Our methodology relies on the prespecification of the distribution F. We will also show that the quantile crossing problem can be overcome by this modelling framework. Equations (2)–(3) state that the conditional distribution of the response at a given value of x has a distribution involving a parameter θ and that the transformed quantile of the distribution becomes a linear predictor of the form (5). This can be achieved by defining link functions that connect (3) to (5). The reason for the linear predictors is that generalized linear modelling [8] is a very well-established method for regression modelling. GLMs are estimated by iteratively reweighted least squares (IRLS) and Fisher scoring, and this algorithm is also adopted by VGLMs and VGAMs.The method presented in this paper differs from conventional quantile regression [1] in that we assume F is known whereas the conventional case does not but use an empirical method instead to obtain the quantiles ξτ: the expectation of the check function ρτ(u) results in the property τ=F(ξτ) which defines the τ-quantile (F is the cumulative distribution function (CDF) of F). In this paper we consider the Fs listed in Table 2.


### 1.2. VGLMs and VGAMs

VGLMs/VGAMs provide the engine and overall modelling framework in this work—theVGAM R package described below fits over 150 models and distributions—therefore we only sketch the details here. VGLMs are defined in terms of M linear predictors, η=η1,…,ηMT, as any statistical model for which the conditional density of y given a d–dimensional vector of explanatory variables, x=(x1,x2,…,xd)T has the form(4)Fy∣x;B=hy,η1,…,ηM;x,for some known function h(·), with B=β1β2⋯βM, a d×M matrix of unknown regression coefficients. Ordinarily, x1≡1 for an intercept.In general, theηj of VGLMs may be applied directly to the M parameters, θj, of any distribution, transformed if necessary, as the jth linear predictor(5)ηj=gjθj=ηjx=βjTx=∑k=1dβjkxk,j=1,…,M,where gj is a VGLM–parameter link function, as in Table 1 (see [6] for further choices) and β(j)k is the kth element of βj. Prior to this work the θj were ‘raw’ parameters such as location, scale, and shape parameters; however, in this present work we define them to be quantiles or a very simple function of quantiles.Table 1
Some VGLM/VGAM link functions. The 4th row is implemented inVGAMextra.
Functions Linksgj(θj) Domain ofθj Link names loge() log ⁡ θ j ( 0 , ∞ ) Logarithmic cloglog() log ⁡ - log ⁡ ( 1 - θ j ) ( 0,1 ) Complementary log–log logit() log ⁡ θ j 1 - θ j ( 0,1 ) Logit logffMeanlink() ‡ logit(θj ) − cloglog(θj ) ( 0,1 ) logffMeanlink rhobit() log ⁡ ⁡ 1 + θ j / ( 1 - θ j ) ( - 1,1 ) rhobit ‡This is the VGLM–link for the mean function of the logarithmic distribution.Table 2
New link functions for the quantiles of some 1–parameter distributions. The selectedG function is also shown.
Distribution θ Support ofy Quantile functionξτ FunctionG Quantile link[η(θ;τ)] Exponential λ ( 0 , ∞ ) - 1 λ log ⁡ ( 1 - τ ) log–link log ⁡ log ⁡ ( 1 - τ ) - 1 / θ Benini s ( y 0 , ∞ ) y 0 exp ⁡ - log ⁡ ( 1 - τ ) s log–link log ⁡ y 0 + log ⁡ ( 1 - τ ) - 1 / θ Rayleigh b ( 0 , ∞ ) b - 2 log ⁡ ( 1 - τ ) log–link log ⁡ θ + 1 2 log ⁡ log ⁡ ( 1 - τ ) - 2 Gamma s ( 0 , ∞ ) No closed–form log–link log ⁡ q g a m m a ( τ , s h a p e = θ ) Maxwell† a ( 0 , ∞ ) 2 a ⋅ q g a m m a ( τ , 1.5 ) ) log–link 1 2 log ⁡ 2  q g a m m a ( τ , 1.5 ) θ Topp–Leone∥ s ( 0,1 ) 1 - 1 - τ 1 / s logit logit 1 - 1 - τ 1 / θ 1–par Normal‡ σ R μ 0 ± 2 ⋅ θ 2 κ ( τ ) ) identity μ 0 ± 2 ⋅ θ 2 κ ( τ ) † qgamma() is the quantile function of the standard gamma distribution in R. ∥ logit ( θ ) = log ⁡ θ / 1 - θ. ‡ κ ( τ ) = e r f - 1 ( 2 τ - 1 ), with erf() denoting the error function.In matrix form one can writeη=η(x)=(6)η1x⋮ηMx=β1Tx⋮βMTx=β11⋯β1p⋮⋱⋮βM1⋯βMpx=∑k=1dβkxk=BTx,where β(k)=β(1)k,β(2)k,…,β(M)kT, k=1,…,d. Sometimes, for some j, it may be required to model ηj asintercept–only, that is, ηj=β(j)1, and β(j)k≡0 for k=2,…,d.VGAMs are a nonparametric extension of VGLMs, that is, (6) is generalized to(7)ηx=β1+∑k=2dfkxk=H1β1∗+∑k=2dHkfk∗xkwith fk∗(xk)=f(1)k∗(xk),…,f(Rk)k∗(xk)T. Usually the component functions are estimated by splines. Here, H1,…,Hd are known full–column rank constraint matrices, and β(1) is a vector of unknown intercepts. With no constraints at all, H1=⋯=Hd=IM (the order-M identity matrix). For VGLMs, the fk are linear so that, cf. (6),(8)B=H1β1⋆∣H2β2⋆∣⋯∣Hdβd⋆.The Hk can enforce a wide range of linear constraints such as parallelism and exchangeability.


### 1.3. Estimation

VGLMs are estimated by maximum likelihood performed by IRLS using the expected information. The VGLM log–likelihood is given by(9)lη=∑i=1nwiliη1xi,…,ηMxi,for known fixed positive prior weights wi, and a Newton–like algorithm for maximizing (9) has the form β(a)=β(a-1)+Iβ(a-1)-1Uβ(a-1), where I is the overall expected information matrix (EIM), U is the score vector, and a is the iteration number. The vector β(a) is obtained as the solution of the generalized least squares problem β(a)=argminβRSS, where the quantity minimized at each IRLS iteration is the weighted (or residual) sum of squares, RSS =(10)∑i=1nwizia-1-ηia-1TWia-1zia-1-ηia-1.The (M×M) Wi are known as the working weight matrices and they have (j,k)th element given by(11)Wij,k=-wiE∂2li∂ηj∂ηk.The use of individual EIMs instead of observed information matrices means that Fisher scoring is used rather than the Newton–Raphson algorithm.VGAMs are also estimated by IRLS, where the difference with respect to VGLMs is that a vector additive model is now fitted to the pseudo–responsezi with explanatory variables xi and working weight matrices Wi at each IRLS iteration. Two approaches are currently used byVGAM to estimate the component functions f∗: regression splines and vector smoothing methods with vector backfitting. Rudimentary P-splines [9] are almost operational, albeit this work is not yet complete. Compared to VGLMs, the VGAM log–likelihood includes a penalty if used with vector smoothing splines. InVGAM the objective function maximized is(12)liη1xi,…,ηMxi-12∑k=1d∑j=1ncolHkλjk∫akbkfjk∗′′t2dt.Here, the λ(j)k are nonnegative smoothing parameters, and ak≤xik≤bk are endpoints covering the values of each covariate. The basic penalty approach adopted here is described in Green and Silverman Green and Silverman [10].


## 1.1. Background

Much of modern regression analysis for estimating conditional quantile functions may be viewed as starting from Koenker and Bassett [1], who offered a systematic strategy for examining how covariates influence the entire response distribution. The fundamental idea is based on the linear specification of the τth quantile function Qy(τ∣x)=βτTx and finding βτ∈Rp that solves the optimization problem(1)minβτ∈Rp⁡∑ρτyi-βτTxi,for independent and identically distributed (i.i.d.) observations from a family of linear quantile regression models, say yi=βτTx+εi,τ,i=1,…,n. Equation (1) can be reformulated as a linear programming problem using the piecewise linear function ρτ(u)=u·τ-I(u<0) for τ∈(0,1). More details can be found in Koenker [2].In the spirit of quantile regression, the conditional distributionY∣x is usually unspecified, although it relies on normal–based asymptotic theory that is used for inference, whilst the assumption of homoskedasticity of the error terms εi,τ is dropped. In this paper we use an alternative approach of conditional–quantile regression based on assuming a prespecified distribution for the response. Parametric quantile regression has some advantages over many nonparametric approaches, including overcoming the quantile crossing problem. Two examples are Noufaily and Jones [3] which is based on the generalized gamma distribution and generalized additive models for location, scale, and shape (GAMLSS; [4]). Further examples are the LMS-BCN method involving the standard normal distribution and a three–parameter Box–Cox transformation [5] and the classical method of quantile regression based on the asymmetric Laplace distribution (ALD).Our approach uses the vector generalized linear and additive model (VGLM/VGAM; [6, 7]) framework. We develop new link functions, G, for the quantile regression model(2)Y∣x~Fx;θ,(3)ητθ=GQyτ∣x,θ,for a vector of quantiles τ=τ1,…,τLT. Our methodology relies on the prespecification of the distribution F. We will also show that the quantile crossing problem can be overcome by this modelling framework. Equations (2)–(3) state that the conditional distribution of the response at a given value of x has a distribution involving a parameter θ and that the transformed quantile of the distribution becomes a linear predictor of the form (5). This can be achieved by defining link functions that connect (3) to (5). The reason for the linear predictors is that generalized linear modelling [8] is a very well-established method for regression modelling. GLMs are estimated by iteratively reweighted least squares (IRLS) and Fisher scoring, and this algorithm is also adopted by VGLMs and VGAMs.The method presented in this paper differs from conventional quantile regression [1] in that we assume F is known whereas the conventional case does not but use an empirical method instead to obtain the quantiles ξτ: the expectation of the check function ρτ(u) results in the property τ=F(ξτ) which defines the τ-quantile (F is the cumulative distribution function (CDF) of F). In this paper we consider the Fs listed in Table 2.


## 1.2. VGLMs and VGAMs

VGLMs/VGAMs provide the engine and overall modelling framework in this work—theVGAM R package described below fits over 150 models and distributions—therefore we only sketch the details here. VGLMs are defined in terms of M linear predictors, η=η1,…,ηMT, as any statistical model for which the conditional density of y given a d–dimensional vector of explanatory variables, x=(x1,x2,…,xd)T has the form(4)Fy∣x;B=hy,η1,…,ηM;x,for some known function h(·), with B=β1β2⋯βM, a d×M matrix of unknown regression coefficients. Ordinarily, x1≡1 for an intercept.In general, theηj of VGLMs may be applied directly to the M parameters, θj, of any distribution, transformed if necessary, as the jth linear predictor(5)ηj=gjθj=ηjx=βjTx=∑k=1dβjkxk,j=1,…,M,where gj is a VGLM–parameter link function, as in Table 1 (see [6] for further choices) and β(j)k is the kth element of βj. Prior to this work the θj were ‘raw’ parameters such as location, scale, and shape parameters; however, in this present work we define them to be quantiles or a very simple function of quantiles.Table 1
Some VGLM/VGAM link functions. The 4th row is implemented inVGAMextra.
Functions Linksgj(θj) Domain ofθj Link names loge() log ⁡ θ j ( 0 , ∞ ) Logarithmic cloglog() log ⁡ - log ⁡ ( 1 - θ j ) ( 0,1 ) Complementary log–log logit() log ⁡ θ j 1 - θ j ( 0,1 ) Logit logffMeanlink() ‡ logit(θj ) − cloglog(θj ) ( 0,1 ) logffMeanlink rhobit() log ⁡ ⁡ 1 + θ j / ( 1 - θ j ) ( - 1,1 ) rhobit ‡This is the VGLM–link for the mean function of the logarithmic distribution.Table 2
New link functions for the quantiles of some 1–parameter distributions. The selectedG function is also shown.
Distribution θ Support ofy Quantile functionξτ FunctionG Quantile link[η(θ;τ)] Exponential λ ( 0 , ∞ ) - 1 λ log ⁡ ( 1 - τ ) log–link log ⁡ log ⁡ ( 1 - τ ) - 1 / θ Benini s ( y 0 , ∞ ) y 0 exp ⁡ - log ⁡ ( 1 - τ ) s log–link log ⁡ y 0 + log ⁡ ( 1 - τ ) - 1 / θ Rayleigh b ( 0 , ∞ ) b - 2 log ⁡ ( 1 - τ ) log–link log ⁡ θ + 1 2 log ⁡ log ⁡ ( 1 - τ ) - 2 Gamma s ( 0 , ∞ ) No closed–form log–link log ⁡ q g a m m a ( τ , s h a p e = θ ) Maxwell† a ( 0 , ∞ ) 2 a ⋅ q g a m m a ( τ , 1.5 ) ) log–link 1 2 log ⁡ 2  q g a m m a ( τ , 1.5 ) θ Topp–Leone∥ s ( 0,1 ) 1 - 1 - τ 1 / s logit logit 1 - 1 - τ 1 / θ 1–par Normal‡ σ R μ 0 ± 2 ⋅ θ 2 κ ( τ ) ) identity μ 0 ± 2 ⋅ θ 2 κ ( τ ) † qgamma() is the quantile function of the standard gamma distribution in R. ∥ logit ( θ ) = log ⁡ θ / 1 - θ. ‡ κ ( τ ) = e r f - 1 ( 2 τ - 1 ), with erf() denoting the error function.In matrix form one can writeη=η(x)=(6)η1x⋮ηMx=β1Tx⋮βMTx=β11⋯β1p⋮⋱⋮βM1⋯βMpx=∑k=1dβkxk=BTx,where β(k)=β(1)k,β(2)k,…,β(M)kT, k=1,…,d. Sometimes, for some j, it may be required to model ηj asintercept–only, that is, ηj=β(j)1, and β(j)k≡0 for k=2,…,d.VGAMs are a nonparametric extension of VGLMs, that is, (6) is generalized to(7)ηx=β1+∑k=2dfkxk=H1β1∗+∑k=2dHkfk∗xkwith fk∗(xk)=f(1)k∗(xk),…,f(Rk)k∗(xk)T. Usually the component functions are estimated by splines. Here, H1,…,Hd are known full–column rank constraint matrices, and β(1) is a vector of unknown intercepts. With no constraints at all, H1=⋯=Hd=IM (the order-M identity matrix). For VGLMs, the fk are linear so that, cf. (6),(8)B=H1β1⋆∣H2β2⋆∣⋯∣Hdβd⋆.The Hk can enforce a wide range of linear constraints such as parallelism and exchangeability.


## 1.3. Estimation

VGLMs are estimated by maximum likelihood performed by IRLS using the expected information. The VGLM log–likelihood is given by(9)lη=∑i=1nwiliη1xi,…,ηMxi,for known fixed positive prior weights wi, and a Newton–like algorithm for maximizing (9) has the form β(a)=β(a-1)+Iβ(a-1)-1Uβ(a-1), where I is the overall expected information matrix (EIM), U is the score vector, and a is the iteration number. The vector β(a) is obtained as the solution of the generalized least squares problem β(a)=argminβRSS, where the quantity minimized at each IRLS iteration is the weighted (or residual) sum of squares, RSS =(10)∑i=1nwizia-1-ηia-1TWia-1zia-1-ηia-1.The (M×M) Wi are known as the working weight matrices and they have (j,k)th element given by(11)Wij,k=-wiE∂2li∂ηj∂ηk.The use of individual EIMs instead of observed information matrices means that Fisher scoring is used rather than the Newton–Raphson algorithm.VGAMs are also estimated by IRLS, where the difference with respect to VGLMs is that a vector additive model is now fitted to the pseudo–responsezi with explanatory variables xi and working weight matrices Wi at each IRLS iteration. Two approaches are currently used byVGAM to estimate the component functions f∗: regression splines and vector smoothing methods with vector backfitting. Rudimentary P-splines [9] are almost operational, albeit this work is not yet complete. Compared to VGLMs, the VGAM log–likelihood includes a penalty if used with vector smoothing splines. InVGAM the objective function maximized is(12)liη1xi,…,ηMxi-12∑k=1d∑j=1ncolHkλjk∫akbkfjk∗′′t2dt.Here, the λ(j)k are nonnegative smoothing parameters, and ak≤xik≤bk are endpoints covering the values of each covariate. The basic penalty approach adopted here is described in Green and Silverman Green and Silverman [10].


## 2. Methodology

LetF(η;y,x) be a 1-parameter statistical model as in (4) parametrized by θ∈Θ⊂R for some parameter space Θ residing in (-∞,∞). Also let Qyτ∣x be the corresponding quantile function with τ∈(0,1). Crucially, note that (5) handles suitable transformations of θ in the linear predictor by parameter link functions. In contrast our proposal focusses on directly modelling Qyτ∣x,θ via a smooth and one–to–one function G, in the form of(13)ητ=GQyτ∣x,θ=G∗θ∣x,τ,which is to be incorporated in the VGLM/VGAM log–likelihood, namely, (9) and (12). Here, τ=τ1,…,τLT is a prespecified vector of quantiles of interest. Examples of (13) are log⁡ξτ, logitξτ, and ξτ.Equation (13) is central to this work. It allows modelling choices via G for the quantile function Qy, and it represents a new modification to the VGLM/VGAM framework. Note that G resembles alink function within the VGLM/VGAM framework as in Table 1. Two notes: first, without any loss of generality, (13) can be seen (strictly) as a function of θ since the quantiles τ and the covariates x are known. Secondly, G∗ is monotonic and one–to–one, as a result of the composite of G and Qy which also hold such properties. However, during the fitting process, the IRLS algorithm internally requires the inverse of G∗. Working with 1–parameter distribution at this stage eases the implementation via Fisher scoring because the inverse G∗-1 can be derived manually and then incorporated in the IRLS algorithm. In a few cases, the inverse of G∗ does not have a closed form, such as the 1–parameter gamma distribution, and an alternative iterative method is employed to approximate G∗-1. To achieve this efficiently, two choices are available. These are (a) newtonRaphson.basic() fromVGAMextra and (b) VGAM::bisection.basic(), twovectorized implementations of the well-known Newton–Raphson and bisection algorithms, to solve the roots of a real-valued function in a given interval (a,b). Further details are given in Section 2.2.One advantage of this work is that the VGLM/VGAM framework can circumvent the quantile crossing problem (e.g., [2, 11], Sect. 2.5) by choosing H1=IM and H2=H3=⋯=Hd=1M (an M-vector of ones). Under this parallelism assumption the method borrows strength across the entire data set so that the additive models for ηj with respect to x2,…,xd are parallel. Each family function in Tables 2 and 3 has a parallel argument which is FALSE by default. Using the syntax ofVGAM based on Chambers and Hastie [12], setting parallel = TRUE (or parallel = FALSE ~ 1) results in H1=IM and H2=H3=⋯=1M; i.e., it is false for only the intercept.Table 3
Inverse of the quantile links and names inVGAMextra. “Approximate” means that Newton–Raphson or bisection is used to approximate the inverse. All family functions except for normal1sdff(), which is in VGAMextra, are in VGAM.
Distribution θ InverseL[θ(η;τ)] Family function Link inVGAMextra Exponential λ - log ⁡ ( 1 - τ ) e η exponential() expQlink() Benini s - log ⁡ ( 1 - τ ) ( η - log ⁡ y 0 ) 2 benini1() benini1Qlink() Rayleigh b exp ⁡ ( η ) - 2 log ⁡ ( 1 - τ ) rayleigh() rayleighQlink() Gamma s Approximate gamma1() gamma1Qlink() Maxwell a 2 q g a m m a τ , 1.5 exp ⁡ ( 2  η ) maxwell() maxwellQlink() Topp–Leone† s log ⁡ τ log ⁡ 1 - [ 1 - l o g i t - 1 ( η ) ] 2 topple() toppleQlink() N(μ=0,σ) σ ( η - μ 0 ) / 2 ⋅ κ ( τ ) normal1sdff() normal1sdQlink() † l o g i t - 1 ( ) denotes the inverse of the logit() transformation.It is noted that for some distributions such as the exponential and Maxwell theηj are naturally parallel with respect to x2,…,xd because log⁡ξτ has the form h1(τ)+h2(x2,…,xd). If this is so then only the intercepts will change as a function of τ and the MLEs for h2 are the same. Other distributions such as the 1-parameter gamma do not possess this property, and then it is necessary to constrain H2=⋯Hd=1M to avoid the quantile crossing problem.


### 2.1. Two Derivations

Ideally the link transforms the support ofY to R because ηj should be unbounded. The three most common cases are as follows. For Y∈(0,∞) a log link is recommended, for Y∈(0,1) a logit link is a good choice, and Y∈(-∞,∞) means that an identity link is natural. These cases have been implemented for seven 1–parameter distributions. The selection of the function G for each Q is shown in the 5th column of Table 2, whilst the resulting quantile links as functions of θ are shown in the last column.We now describe the quantile links for the exponential and the Topp–Leone distributions as examples. Firstly, forY~exponential(θ), with a rate parameter θ>0, the density and CDF are given by f(y;θ)=θe-θy and F(y;θ)=1-e-θy. With a slight change in notation, the quantile function is given by F-1, i.e.,(14)QYτ;θ=-1θlog⁡1-τ,which lies in (0,∞) regardless of the values of τ and θ. Given that values of τ are known (prespecified by the user), (14) becomes a function of θ. Thus, the new quantile link for the exponential distribution as shown in Table 2 is simply obtained by taking G as the logarithmic transformation, as follows:(15)ηθ;τ=log⁡-1θlog⁡1-τ=log⁡log⁡1-τ-1/θ.This quantile link has been implemented inVGAMextra via the function expQlink(), as shown in Table 3. Its inverse (denoted as θ(η;τ)) can be manually obtained from the inverse of (15). Note that the corresponding family function (exponential()) implemented inVGAM includes a (known) location parameter A, which gives the density f(y;θ)=θe-θ(y-A). By default A=0, and it is handled by the argument location.Secondly, consider the Topp–Leone distributionY~Topp-Leone(s) whose support is (0,1) and(16)QYs;τj=1-1-τj1/s,with 0<τj<1. Here, θ=s. To verify this restriction note that 1/s>1, for any shape parameter s∈(0,1), and hence for any τj∈(0,1),(17)0<τj1/s<1⟺0<1-τj1/s<1⟺0<1-1-τj1/s<1.Thus, to allow the quantile function to be modelled by covariates, we take the logit transformation as G. The resulting quantile link for this distribution is simply η(s;τ)=log⁡QY(s;τ), shown in Table 2. The distribution has CDF F(y;s)=y·(2-y)s for 0<y<1, and density f(y;s)=2s(1-y)·[y(2-y)]s-1. The quantile function derives from solving the equation τ0=F(y;s)=y·(2-y)s, for 0<τ0<1, which leads to the quadratic equation y2-2y+τ01/s=0. The solution must lie in (0,1) and is in fact (16), as a function of s. The family function topple() fromVGAM estimates s, where the default link is η(s)=logit(s).


### 2.2. Software Implementation

For practical use by others, we have implemented seven VGLM–quantile links,ητ in theR packageVGAMextra. They are summarized in Table 2. The packageVGAM is a requirement ofVGAMextra because the modelling functions vglm() and vgam(), and all but the last family function of Table 2, reside there. For this paperVGAMextra 0.0-2 andVGAM 1.1-0 or later are required; they are available at www.stat.auckland.ac.nz/~vmir178 and www.stat.auckland.ac.nz/~yee/VGAM/prerelease/ whilst older versions of both are available on CRAN (http://CRAN.R-project.org).One special case isgamma1Qlink(), for the 1–parameter (shape) gamma distribution, defined as(18)ηθ;τ=log⁡gammaτ,shape=θ,whose primary arguments are τ and θ. Its inverse (Table 3) does not admit a closed form and it is approximated by the function VGAM::newtonRaphson.basic(), avectorized implementation of the Newton–Raphson algorithm. Almost all implementations elsewhere of this are for a scalar argument, but we operate on vectors of length n. It works as follows. Our data is effectively yi={yi,xi,d},i=1,…,n, whilst the quantiles of interest, τ or p, must be entered by the user. The shape parameter θ is estimated by IRLS and therefore it is available at each iteration. Thus, for each η(θ;p)0, the ‘inverse’ is given by the root, θ, of the function(19)fθ;p,η=ηθ;p0-log⁡gammaτ,shape=θ.Finally, the inverse of all the VGLM–quantile links is shown in Table3, as well as the name of the corresponding implementation inVGAMextra. The inverse–links are required at different stages of the IRLS by Fisher scoring, which internally switches between η(θ;τ) (namely, Table 2) and θ(η;τ) (namely, Table 3). Specifically, the algorithm requires the score vector and the EIMs at each IRLS iteration, which are given by the following chain–rule formulas:(20)∂l∂η=∂l∂θ·∂θ∂η,-E∂2l∂η2=-E∂2l∂θ2∂θ∂η2.Internally, the functions utilized to compute the inverse areVGAM::eta2theta() or VGAM::theta2eta(). TheVGAMextra Manual and Miranda-Soberanis [13] give further details about the derivation of the quantile links, whilst Yee [6] describes in the IRLS and Fisher scoring algorithms for estimating VGLMs and VGAMs. Complements at the second author’s homepage give additional details on link functions.


### 2.3. Software Use

For the user, this methodology runs as usual by calling the modelling functionsVGAM::vglm() and VGAM::vgam(), except for two modifications that are described below.To start with, we give the following output that shows the central arguments handled byVGAM::vglm():The first adjustment takes place with the argumentformula, a symbolic description of the model to be fit. Usually, an expression like y ~ x2 + x3 should suffice for a response y and covariates x2 and x3. This effectively works for univariate and even for multiple responses say y1, y2, and y3, where the only change is to set cbind(y1, y2, y3)~x2 + x3. Here, the right hand side (RHS) of the formula is applied to each linear predictor.For quantile modelling using VGLMs and VGAMs,Q.reg() must be incorporated in the formula, whose arguments are shown in Table 4. For a given set of quantiles of interest, entered through τ=(τ1,…,τL)T, Q.reg() replicates the response matrixY into NOS·dim(τ) columns, where NOS denotes the number of columns ofY. Then, the RHS of the formula applies to every set of columns according to the number of quantiles of interest. Ordinarily the response is a vector so that NOS=1 and L=M.Table 4
Arguments handled by the functionVGAMextra::Q.reg().
Argument Description y Numeric, a vector or a matrix. It is the response or dependent variable in the formula of the model to be fitted.  pvector A prototype vector. Entries are the conditional p-quantiles in the fitting process.  length.arg A unit-length positive integer. It is the number of p-quantiles to be modelled.As an example, suppose that we have two responsesY1 and Y2 sampled from a prespecified distribution F, as per Table 3, and the quantiles of interest are p=(0.25,0.50,0.75)T. Then Q.reg(cbind(Y1, Y2), pvector = p) will return a matrix with six columns, with the first three columns being Y1, one for each quantile, and similarly the last three columns equal Y2. Thus vglm() handles this model as a multiple responses fit.The second adjustment is related to the argumentfamily, a function that describes the statistical model to be fitted. Each family has at least one argument for the link functions to be used in the fitting process (the name changes from family to family). For example, for VGAM::exponential() this is called link, whilst for the family function VGAM::benini1() (see the third column of Table 3), it is called lshape. When VGLM–quantile modelling is to be performed, the corresponding link (last column of Table 3) must be entered into the family accordingly. All the quantile links manage the same arguments, including p, the vector of quantiles, except by benini1Qlink() which has the additional argument y0.With both modifications, a typical call has the following form:Further fitting variants can be incorporated here, e.g., categorical covariates and the use of smoothers such as regression splines. These and a few other features are illustrated in the following section.


## 2.1. Two Derivations

Ideally the link transforms the support ofY to R because ηj should be unbounded. The three most common cases are as follows. For Y∈(0,∞) a log link is recommended, for Y∈(0,1) a logit link is a good choice, and Y∈(-∞,∞) means that an identity link is natural. These cases have been implemented for seven 1–parameter distributions. The selection of the function G for each Q is shown in the 5th column of Table 2, whilst the resulting quantile links as functions of θ are shown in the last column.We now describe the quantile links for the exponential and the Topp–Leone distributions as examples. Firstly, forY~exponential(θ), with a rate parameter θ>0, the density and CDF are given by f(y;θ)=θe-θy and F(y;θ)=1-e-θy. With a slight change in notation, the quantile function is given by F-1, i.e.,(14)QYτ;θ=-1θlog⁡1-τ,which lies in (0,∞) regardless of the values of τ and θ. Given that values of τ are known (prespecified by the user), (14) becomes a function of θ. Thus, the new quantile link for the exponential distribution as shown in Table 2 is simply obtained by taking G as the logarithmic transformation, as follows:(15)ηθ;τ=log⁡-1θlog⁡1-τ=log⁡log⁡1-τ-1/θ.This quantile link has been implemented inVGAMextra via the function expQlink(), as shown in Table 3. Its inverse (denoted as θ(η;τ)) can be manually obtained from the inverse of (15). Note that the corresponding family function (exponential()) implemented inVGAM includes a (known) location parameter A, which gives the density f(y;θ)=θe-θ(y-A). By default A=0, and it is handled by the argument location.Secondly, consider the Topp–Leone distributionY~Topp-Leone(s) whose support is (0,1) and(16)QYs;τj=1-1-τj1/s,with 0<τj<1. Here, θ=s. To verify this restriction note that 1/s>1, for any shape parameter s∈(0,1), and hence for any τj∈(0,1),(17)0<τj1/s<1⟺0<1-τj1/s<1⟺0<1-1-τj1/s<1.Thus, to allow the quantile function to be modelled by covariates, we take the logit transformation as G. The resulting quantile link for this distribution is simply η(s;τ)=log⁡QY(s;τ), shown in Table 2. The distribution has CDF F(y;s)=y·(2-y)s for 0<y<1, and density f(y;s)=2s(1-y)·[y(2-y)]s-1. The quantile function derives from solving the equation τ0=F(y;s)=y·(2-y)s, for 0<τ0<1, which leads to the quadratic equation y2-2y+τ01/s=0. The solution must lie in (0,1) and is in fact (16), as a function of s. The family function topple() fromVGAM estimates s, where the default link is η(s)=logit(s).


## 2.2. Software Implementation

For practical use by others, we have implemented seven VGLM–quantile links,ητ in theR packageVGAMextra. They are summarized in Table 2. The packageVGAM is a requirement ofVGAMextra because the modelling functions vglm() and vgam(), and all but the last family function of Table 2, reside there. For this paperVGAMextra 0.0-2 andVGAM 1.1-0 or later are required; they are available at www.stat.auckland.ac.nz/~vmir178 and www.stat.auckland.ac.nz/~yee/VGAM/prerelease/ whilst older versions of both are available on CRAN (http://CRAN.R-project.org).One special case isgamma1Qlink(), for the 1–parameter (shape) gamma distribution, defined as(18)ηθ;τ=log⁡gammaτ,shape=θ,whose primary arguments are τ and θ. Its inverse (Table 3) does not admit a closed form and it is approximated by the function VGAM::newtonRaphson.basic(), avectorized implementation of the Newton–Raphson algorithm. Almost all implementations elsewhere of this are for a scalar argument, but we operate on vectors of length n. It works as follows. Our data is effectively yi={yi,xi,d},i=1,…,n, whilst the quantiles of interest, τ or p, must be entered by the user. The shape parameter θ is estimated by IRLS and therefore it is available at each iteration. Thus, for each η(θ;p)0, the ‘inverse’ is given by the root, θ, of the function(19)fθ;p,η=ηθ;p0-log⁡gammaτ,shape=θ.Finally, the inverse of all the VGLM–quantile links is shown in Table3, as well as the name of the corresponding implementation inVGAMextra. The inverse–links are required at different stages of the IRLS by Fisher scoring, which internally switches between η(θ;τ) (namely, Table 2) and θ(η;τ) (namely, Table 3). Specifically, the algorithm requires the score vector and the EIMs at each IRLS iteration, which are given by the following chain–rule formulas:(20)∂l∂η=∂l∂θ·∂θ∂η,-E∂2l∂η2=-E∂2l∂θ2∂θ∂η2.Internally, the functions utilized to compute the inverse areVGAM::eta2theta() or VGAM::theta2eta(). TheVGAMextra Manual and Miranda-Soberanis [13] give further details about the derivation of the quantile links, whilst Yee [6] describes in the IRLS and Fisher scoring algorithms for estimating VGLMs and VGAMs. Complements at the second author’s homepage give additional details on link functions.


## 2.3. Software Use

For the user, this methodology runs as usual by calling the modelling functionsVGAM::vglm() and VGAM::vgam(), except for two modifications that are described below.To start with, we give the following output that shows the central arguments handled byVGAM::vglm():The first adjustment takes place with the argumentformula, a symbolic description of the model to be fit. Usually, an expression like y ~ x2 + x3 should suffice for a response y and covariates x2 and x3. This effectively works for univariate and even for multiple responses say y1, y2, and y3, where the only change is to set cbind(y1, y2, y3)~x2 + x3. Here, the right hand side (RHS) of the formula is applied to each linear predictor.For quantile modelling using VGLMs and VGAMs,Q.reg() must be incorporated in the formula, whose arguments are shown in Table 4. For a given set of quantiles of interest, entered through τ=(τ1,…,τL)T, Q.reg() replicates the response matrixY into NOS·dim(τ) columns, where NOS denotes the number of columns ofY. Then, the RHS of the formula applies to every set of columns according to the number of quantiles of interest. Ordinarily the response is a vector so that NOS=1 and L=M.Table 4
Arguments handled by the functionVGAMextra::Q.reg().
Argument Description y Numeric, a vector or a matrix. It is the response or dependent variable in the formula of the model to be fitted.  pvector A prototype vector. Entries are the conditional p-quantiles in the fitting process.  length.arg A unit-length positive integer. It is the number of p-quantiles to be modelled.As an example, suppose that we have two responsesY1 and Y2 sampled from a prespecified distribution F, as per Table 3, and the quantiles of interest are p=(0.25,0.50,0.75)T. Then Q.reg(cbind(Y1, Y2), pvector = p) will return a matrix with six columns, with the first three columns being Y1, one for each quantile, and similarly the last three columns equal Y2. Thus vglm() handles this model as a multiple responses fit.The second adjustment is related to the argumentfamily, a function that describes the statistical model to be fitted. Each family has at least one argument for the link functions to be used in the fitting process (the name changes from family to family). For example, for VGAM::exponential() this is called link, whilst for the family function VGAM::benini1() (see the third column of Table 3), it is called lshape. When VGLM–quantile modelling is to be performed, the corresponding link (last column of Table 3) must be entered into the family accordingly. All the quantile links manage the same arguments, including p, the vector of quantiles, except by benini1Qlink() which has the additional argument y0.With both modifications, a typical call has the following form:Further fitting variants can be incorporated here, e.g., categorical covariates and the use of smoothers such as regression splines. These and a few other features are illustrated in the following section.


## 3. Examples

### 3.1. Maxwell Data

We use simulation to generaten=200 random variates from a Maxwell distribution whose rate parameter is a function of a single covariate x2. To account for a nonlinear trend in the dataset, additive models with cubic smoothing splines appear to be a better choice over linear schemes such as with VGLMs. In this example we perform the following steps to confirm the performance of the methodology.(1)
Generate random deviates from the Maxwell distribution.(2)
Run conditional VGAM–quantile modelling usingmaxwellQlink() based on the VGAM family function VGAM::maxwell(), which estimates the Maxwell distribution by Fisher scoring.(3)
Perform ordinary quantile regression usingVGAM::alaplace1() that estimates the 1–parameter ALD by Fisher scoring. Here, the special argument tau will be employed.(4)
Plot the artificial data with the estimated quantile functions,Qyτ∣x,θ^ (from (2)), and the estimated quantile curves (from (3)) superimposed.We will consider the quantiles 25%, 50%, and 75% for simplicity, so thatτ=(1/4,1/2,3/4)T.Regarding (1), the data is generated byVGAM::rmaxwell(), which gives random deviates from the Maxwell distribution whose density is f(y;a)=2/πa3/2y2exp(-ay2/2). We use the rate function(21)a=exp⁡2-6sin⁡2xi2-1/5xi2+1/22,where Xi2~i.i.dUnif(0,1), i=1,…,n. The following code chunk sets things up and the dataset is saved as maxdata.The following code chunk performs steps (2) and (3). Note the fitting of additive models viaVGAM::vgam() with smooth terms defined by VGAM::s() where x2 is to be smoothed. To compare both fits, they are saved in fit.Qmodelling (from (2)) and fit.Qregression (from (3)).Figure1 shows the simulated data, the estimated quantile functions, and the fitted quantile curves from fit.Qmodelling and fit.Qregression, obtained from vector smoothing spline fits [14]. The results are similar for x2>0.3, but our present work performs better at the bottom LHS tail. The data coverage from each modelling framework is summarized in Table 5. Once again our work outperforms the ALD method.Table 5
Empirical data coverage from quantile–modelling usingVGAMextra (QM–VGAMextra) and quantile–regression from VGAM (QR–VGAM), after fitting (21).
τj QM–VGAMextra coverage QR–VGAM coverage 25% 26% 28.5% 50% 50% 54% 75% 73.5% 78.5%Figure 1
Simulated Maxwell data (21) including (a) the fitted quantile functions from fit.Qmodelling (VGAMextra) and (b) the fitted quantile curves from fit.Qregression (VGAM). The quantile curves in both cases derive from vector smoothing spline fits.We conclude with a few remarks.(1)
The argumentp is available forall quantile links in Table 3 and not only for maxwellQlink(). It can be assigned any vector of percentile values.(2)
Under the conditional VGAM–quantile modelling framework, the arguments to handle the parallelism assumption such as the argumentsparallel.locat and parallel.scale in family functions are no longer required. This is internally managed by the new quantile links rather than being managed by the family function.(3)
Iffit is a Qlink fit, then fitted(fit) returns the fitted quantiles. This is in the form of a n×(L·NOS) matrix. Similarly, predict(fit) returns a n×M matrix where the ith row is ηiT.


### 3.2. Comparison with the Quantreg Package

For checking purposes, the results are compared withquantreg too. Figure 2 gives the results based on the following code.Figure 2
Simulated Maxwell data (21) including (a) the fitted quantile functions from fit.Qmodelling (VGAMextra) and (b) the fitted quantile curves from qrfit50, etc. fromquantreg.The results should be the similar to Section3.1 because the ALD and the classical quantile regression method are essentially the same. It can be seen that the bottom LHS corner is not modelled well withquantreg either. Once again our method performs best, which is not surprising given the strong distributional assumption.


### 3.3. Exponential Data

Feigl and Zelen [15] fit an exponential distribution to a data set comprising the time to death (in weeks) and white blood cell counts for two groups of leukaemia patients, and a binary variable for AG-positive and AG-negative. The two groups were not created by random allocation. The variable AG is the morphological variable, the AG factor; a numeric vector where 1 means AG-positive and 2 means AG-negative. We create AG01 which is AG - 1. We take the log of the white blood cell count (WBC) because it is very highly skewed. The data are found inGLMsData on CRAN, which supports Dunn and Smyth [16].One benefit of quantile modelling with VGLMs is that it easily allows comparisons of the effect ofAG01 or any other indicator variable, at different quantiles. First note that, for AG-positive patients with logWBC = 9, the 25% percentile for time to death is around exp⁡2.848≈17.25 weeks, whilst the 75% percentile is about exp⁡4.42≈83.1 weeks. Secondly, the coefficient of AG011 measures the influence of the AG factor on the time to death. Keeping the levels of WBC constant, for patients at either the 25% or 75% percentiles, the time to death for AG–negatives compared to AG–positives is multiplicative by a factor of exp⁡(-1.02)≈0.361, i.e., a 63.9% reduction in lifetime.For further illustration’s sake, we fit a 1-parameter gamma distribution to these data and interpret the results. Unlike the Maxwell and exponential distributions, where simple mathematics shows that different quantiles are parallel because their logarithm is additive with respect toτ, the 1-parameter gamma does not possess this property.Here, keeping the level of WBC constant, for patients at the 25% percentile, the time to death for AG–negatives compared to AG–positives is multiplicative by a factor ofexp⁡(-1.46)≈0.232, i.e., a 76.8% reduction. In comparison, for patients at the 75% percentile, the time to death for AG–negatives compared to AG–positives is multiplicative by a factor of exp⁡(-1.27)≈0.281, i.e., a 71.9% reduction. This suggests that the effect of AG is greater for more severe cases than those who live longer in general.Finally, just to check, we obtain the constraint matrices for each predictor:There is a parallelism assumption made forlogWBC but not for any of the other explanatory variables.


## 3.1. Maxwell Data

We use simulation to generaten=200 random variates from a Maxwell distribution whose rate parameter is a function of a single covariate x2. To account for a nonlinear trend in the dataset, additive models with cubic smoothing splines appear to be a better choice over linear schemes such as with VGLMs. In this example we perform the following steps to confirm the performance of the methodology.(1)
Generate random deviates from the Maxwell distribution.(2)
Run conditional VGAM–quantile modelling usingmaxwellQlink() based on the VGAM family function VGAM::maxwell(), which estimates the Maxwell distribution by Fisher scoring.(3)
Perform ordinary quantile regression usingVGAM::alaplace1() that estimates the 1–parameter ALD by Fisher scoring. Here, the special argument tau will be employed.(4)
Plot the artificial data with the estimated quantile functions,Qyτ∣x,θ^ (from (2)), and the estimated quantile curves (from (3)) superimposed.We will consider the quantiles 25%, 50%, and 75% for simplicity, so thatτ=(1/4,1/2,3/4)T.Regarding (1), the data is generated byVGAM::rmaxwell(), which gives random deviates from the Maxwell distribution whose density is f(y;a)=2/πa3/2y2exp(-ay2/2). We use the rate function(21)a=exp⁡2-6sin⁡2xi2-1/5xi2+1/22,where Xi2~i.i.dUnif(0,1), i=1,…,n. The following code chunk sets things up and the dataset is saved as maxdata.The following code chunk performs steps (2) and (3). Note the fitting of additive models viaVGAM::vgam() with smooth terms defined by VGAM::s() where x2 is to be smoothed. To compare both fits, they are saved in fit.Qmodelling (from (2)) and fit.Qregression (from (3)).Figure1 shows the simulated data, the estimated quantile functions, and the fitted quantile curves from fit.Qmodelling and fit.Qregression, obtained from vector smoothing spline fits [14]. The results are similar for x2>0.3, but our present work performs better at the bottom LHS tail. The data coverage from each modelling framework is summarized in Table 5. Once again our work outperforms the ALD method.Table 5
Empirical data coverage from quantile–modelling usingVGAMextra (QM–VGAMextra) and quantile–regression from VGAM (QR–VGAM), after fitting (21).
τj QM–VGAMextra coverage QR–VGAM coverage 25% 26% 28.5% 50% 50% 54% 75% 73.5% 78.5%Figure 1
Simulated Maxwell data (21) including (a) the fitted quantile functions from fit.Qmodelling (VGAMextra) and (b) the fitted quantile curves from fit.Qregression (VGAM). The quantile curves in both cases derive from vector smoothing spline fits.We conclude with a few remarks.(1)
The argumentp is available forall quantile links in Table 3 and not only for maxwellQlink(). It can be assigned any vector of percentile values.(2)
Under the conditional VGAM–quantile modelling framework, the arguments to handle the parallelism assumption such as the argumentsparallel.locat and parallel.scale in family functions are no longer required. This is internally managed by the new quantile links rather than being managed by the family function.(3)
Iffit is a Qlink fit, then fitted(fit) returns the fitted quantiles. This is in the form of a n×(L·NOS) matrix. Similarly, predict(fit) returns a n×M matrix where the ith row is ηiT.


## 3.2. Comparison with the Quantreg Package

For checking purposes, the results are compared withquantreg too. Figure 2 gives the results based on the following code.Figure 2
Simulated Maxwell data (21) including (a) the fitted quantile functions from fit.Qmodelling (VGAMextra) and (b) the fitted quantile curves from qrfit50, etc. fromquantreg.The results should be the similar to Section3.1 because the ALD and the classical quantile regression method are essentially the same. It can be seen that the bottom LHS corner is not modelled well withquantreg either. Once again our method performs best, which is not surprising given the strong distributional assumption.


## 3.3. Exponential Data

Feigl and Zelen [15] fit an exponential distribution to a data set comprising the time to death (in weeks) and white blood cell counts for two groups of leukaemia patients, and a binary variable for AG-positive and AG-negative. The two groups were not created by random allocation. The variable AG is the morphological variable, the AG factor; a numeric vector where 1 means AG-positive and 2 means AG-negative. We create AG01 which is AG - 1. We take the log of the white blood cell count (WBC) because it is very highly skewed. The data are found inGLMsData on CRAN, which supports Dunn and Smyth [16].One benefit of quantile modelling with VGLMs is that it easily allows comparisons of the effect ofAG01 or any other indicator variable, at different quantiles. First note that, for AG-positive patients with logWBC = 9, the 25% percentile for time to death is around exp⁡2.848≈17.25 weeks, whilst the 75% percentile is about exp⁡4.42≈83.1 weeks. Secondly, the coefficient of AG011 measures the influence of the AG factor on the time to death. Keeping the levels of WBC constant, for patients at either the 25% or 75% percentiles, the time to death for AG–negatives compared to AG–positives is multiplicative by a factor of exp⁡(-1.02)≈0.361, i.e., a 63.9% reduction in lifetime.For further illustration’s sake, we fit a 1-parameter gamma distribution to these data and interpret the results. Unlike the Maxwell and exponential distributions, where simple mathematics shows that different quantiles are parallel because their logarithm is additive with respect toτ, the 1-parameter gamma does not possess this property.Here, keeping the level of WBC constant, for patients at the 25% percentile, the time to death for AG–negatives compared to AG–positives is multiplicative by a factor ofexp⁡(-1.46)≈0.232, i.e., a 76.8% reduction. In comparison, for patients at the 75% percentile, the time to death for AG–negatives compared to AG–positives is multiplicative by a factor of exp⁡(-1.27)≈0.281, i.e., a 71.9% reduction. This suggests that the effect of AG is greater for more severe cases than those who live longer in general.Finally, just to check, we obtain the constraint matrices for each predictor:There is a parallelism assumption made forlogWBC but not for any of the other explanatory variables.


## 4. Discussion and Future Work

This work in parametric quantile regression is blighted by the strong assumption of the assumed distribution. In theory, this might be ameliorated somewhat by implementing as many distributions as possible. Some of the distributions listed in Table2 have real applications, for example, in kinetic-molecular theory the speed of individual molecules of idealized gases follows the Maxwell distribution and the average kinetic speed is directly related to Kelvin temperature. In experiments that do not satisfy the various postulates made (such as the effects of the container) one might model the median particle speed with x comprising temperature and other covariates such as volume of the container and density of the container walls. Different forms of gases, such as plasmas and rarefied gases, could be modelled as such too. Another example is the Rayleigh distribution which is similar to the Maxwell distribution. In two-dimensions, and in applications of magnetic resonance imaging (MRI), complex images are often viewed in terms of the background data, which is Rayleigh distributed. Nonstandard background information could be included in x and their effects on the distribution examined.In the current software implementation there are limitations due to its internal design. For example, it would be good ifworked like many otherVGAM models. The difficulty here is that the @linkinv S4 slot of aVGAM family function has eta as an argument, and in our implementation this could only possibly be created by supplying the new percentiles to predict() beforehand.Another minor deficiency in our software implementation is that the response vector is replicateddim⁡(τ) times so that is a form of recycling. Possibly this could be avoided because the memory requirement might be excessive when either dim⁡(τ) or n are very large.At present, the VGAM framework has infrastructure to afford 1–parameter quantile links. For quantile functions depending on 2 or more parameter, such as the two-parameter gamma distribution, the quantiles will be bivariate functions whose inverse would probably not admit a closed form. Nevertheless, future work includes being able to write links for two-parameter distributions, of which the normal distribution would be the most important. For this, the methodology behind Yee and Miranda-Soberanis [17] could be employed; they solve a decades-old problem implementing the two-parameter canonical link function log⁡(μ/(μ+k)) of the negative binomial distribution. We have already commenced work in this direction, e.g., with the 2-parameter gamma distribution.



---

*Source: 3493628-2019-05-07.xml*