Package 'quantreg'

Description

Univariate adaptive kernel density estimation a la Silverman. As used by Portnoy and Koenker (1989).

Usage

akj(x, z =, p =, h = -1, alpha = 0.5, kappa = 0.9, iker1 = 0)

Arguments

Value

a list structure is with components

Note

if the score function values are of interest, the Cauchy kernel may be preferable.

References

Portnoy, S and R Koenker, (1989) Adaptive L Estimation of Linear Models; Annals of Statistics 17, 362–81.

Silverman, B. (1986) Density Estimation, pp 100–104.

Examples

set.seed(1)
 x <- c(rnorm(600), 2 + 2*rnorm(400))
 xx <- seq(-5, 8, length=200)
 z <- akj(x, xx)
 plot(xx, z$dens, ylim=range(0,z$dens), type ="l", col=2)
 abline(h=0, col="gray", lty=3)
 plot(xx, z$psi, type ="l", col=2, main = expression(hat(psi(x))))
 plot(xx, z$score, type ="l", col=2,
      main = expression("score " * hat(psi) * "'" * (x)))

 if(require("nor1mix")) {
  m3 <- norMix(mu= c(-4, 0, 3), sigma = c(1/3, 1, 2),
               w = c(.1,.5,.4))
  plot(m3, p.norm = FALSE)
  set.seed(11)
  x <- rnorMix(1000, m3)
  z2 <- akj(x, xx)
  lines(xx, z2$dens, col=2)
  z3 <- akj(x, xx, kappa = 0.5, alpha = 0.88)
  lines(xx, z3$dens, col=3)
 }

Description

Compute test statistics for two or more quantile regression fits.

Usage

## S3 method for class 'rq'
anova(object, ..., test = "Wald", joint = TRUE, score =
                       "tau", se = "nid", iid = TRUE, R = 200, trim = NULL)
## S3 method for class 'rqs'
anova(object, ..., se = "nid", iid = TRUE, joint = TRUE)
## S3 method for class 'rqlist'
anova(object, ...,  test = "Wald", joint = TRUE, 
	score = "tau", se = "nid", iid = TRUE, R = 200, trim = NULL)
rq.test.rank(x0, x1, y, v = NULL, score = "wilcoxon", weights = NULL, tau=.5,
        iid = TRUE, delta0 = rep(0,NCOL(x1)), omega = 1, trim = NULL, pvalue = "F")
rq.test.anowar(x0,x1,y,tau,R)
## S3 method for class 'anova.rq'
print(x, ...)

Arguments

Details

There are two (as yet) distinct forms of the test. In the first the fitted objects all have the same specified quantile (tau) and the intent is to test the hypothesis that smaller models are adequate relative to the largest specified model. In the second form of the test the linear predictor of the fits are all the same, but the specified quantiles (taus) are different.

In the former case there are three options for the argument ‘test’, by default a Wald test is computed as in Bassett and Koenker (1982). If test = 'anowar' is specified then the test is based on the procedure suggested in Chen, Ying, Zhang and Zhao (2008); the test is based on the difference in the QR objective functions at the restricted and unrestricted models with a reference distribution computed by simulation. The p-value of this form of the test is produced by fitting a density to the simulation values forming the reference distribution using the logspline function from the logspline package. The acronym anowar stands for analysis of weighted absolute residuals. If test='rank' is specified, then a rank test statistic is computed as described in Gutenbrunner, Jureckova, Koenker and Portnoy (1993). In the latter case one can also specify a form for the score function of the rank test, by default the Wilcoxon score is used, the other options are score=‘sign’ for median (sign) scores, or score=‘normal’ for normal (van der Waerden) scores. A fourth option is score=‘tau’ which is a generalization of median scores to an arbitrary quantile, in this case the quantile is assumed to be the one associated with the fitting of the specified objects. The computing of the rank form of the test is carried out in the rq.test.rank function, see ranks for further details on the score function options. The Wald form of the test is local in sense that the null hypothesis asserts only that a subset of the covariates are “insignificant” at the specified quantile of interest. The rank form of the test can also be used to test the global hypothesis that a subset is “insignificant” over an entire range of quantiles. The use of the score function score = "tau" restricts the rank test to the local hypothesis of the Wald test.

In the latter case the hypothesis of interest is that the slope coefficients of the models are identical. The test statistic is a variant of the Wald test described in Koenker and Bassett (1982).

By default, both forms of the tests return an F-like statistic in the sense that the an asymptotically Chi-squared statistic is divided by its degrees of freedom and the reported p-value is computed for an F statistic based on the numerator degrees of freedom equal to the rank of the null hypothesis and the denominator degrees of freedom is taken to be the sample size minus the number of parameters of the maintained model.

Value

An object of class ‘"anova"’ inheriting from class ‘"data.frame"’.

WARNING

An attempt to verify that the models are nested in the first form of the test is made, but this relies on checking set inclusion of the list of variable names and is subject to obvious ambiguities when variable names are generic. The comparison between two or more models will only be valid if they are fitted to the same dataset. This may be a problem if there are missing values and R's default of ‘na.action = na.omit’ is used. The rank version of the nested model tests involves computing the entire regression quantile process using parametric linear programming and thus can be rather slow and memory intensive on problems with more than several thousand observations.

Author(s)

Roger Koenker

References

[1] Bassett, G. and R. Koenker (1982). Tests of Linear Hypotheses and L1 Estimation, Econometrica, 50, 1577–83.

[2] Koenker, R. W. and Bassett, G. W. (1982). Robust Tests for Heteroscedasticity based on Regression Quantiles, Econometrica, 50, 43–61.

[3] Gutenbrunner, C., Jureckova, J., Koenker, R, and S. Portnoy (1993). Tests of Linear Hypotheses based on Regression Rank Scores, J. of Nonparametric Statistics, 2, 307–331.

[4] Chen, K. Z. Ying, H. Zhang, and L Zhao, (2008) Analysis of least absolute deviations, Biometrika, 95, 107-122.

[5] Koenker, R. W. (2005). Quantile Regression, Cambridge U. Press.

See Also

The model fitting function rq, and the functions for testing hypothesis on the entire quantile regression process KhmaladzeTest. For further details on the rank tests see ranks.

Examples

data(barro)
fit0 <- rq(y.net ~  lgdp2 + fse2 + gedy2 , data = barro)
fit1 <- rq(y.net ~  lgdp2 + fse2 + gedy2 + Iy2 + gcony2, data = barro)
fit2 <- rq(y.net ~  lgdp2 + fse2 + gedy2 + Iy2 + gcony2, data = barro,tau=.75)
fit3 <- rq(y.net ~  lgdp2 + fse2 + gedy2 + Iy2 + gcony2, data = barro,tau=.25)
anova(fit1,fit0)
anova(fit1,fit2,fit3)
anova(fit1,fit2,fit3,joint=FALSE)
# Alternatively fitting can be done in one call:
fit <- rq(y.net ~  lgdp2 + fse2 + gedy2 + Iy2 + gcony2, 
	  method = "fn", tau = 1:4/5, data = barro)

Description

function to compute bandwidth for sparsity estimation

Usage

bandwidth.rq(p, n, hs=TRUE, alpha=0.05)

Arguments

Details

If hs=TRUE (default) then the Hall-Sheather(1988) rule $O(n^{-1/3})$ is used, if hs=FALSE then the Bofinger $O(n^{-1/5})$ is used.

Value

returns a vector of bandwidths corresponding to the argument p.

Author(s)

Roger Koenker [email protected]

References

Hall and Sheather(1988, JRSS(B)),Bofinger (1975, Aus. J. Stat)

Description

Version of the Barro Growth Data used in Koenker and Machado(1999). This is a regression data set consisting of 161 observations on determinants of cross country GDP growth rates. There are 13 covariates with dimnames corresponding to the original Barro and Lee source. See https://www.nber.org/pub/barro.lee/. The first 71 observations are on the period 1965-75, remainder on 1987-85.

Usage

data(barro)

Format

A data frame containing 161 observations on 14 variables:

References

Koenker, R. and J.A.F. Machado (1999) Goodness of Fit and Related Inference Processes for Quantile Regression, JASA, 1296-1310.

Description

Functions used to estimated standard errors, confidence intervals and tests of hypotheses for censored quantile regression models using the Portnoy and Peng-Huang methods.

Usage

boot.crq(x, y, c, taus, method, ctype = "right", R = 100, mboot, bmethod = "jack", ...)

Arguments

Details

There are several refinements that are still unimplemented. Percentile methods should be incorporated, and extensions of the methods to be used in anova.rq should be made. Note that bootstrapping for the Powell method "Powell" is done via boot.rq. For problems with n > 3000 a message is printed indicated progress in the resampling.

Value

A matrix of dimension R by p is returned with the R resampled estimates of the vector of quantile regression parameters. When mofn < n for the "xy" method this matrix has been deflated by the factor sqrt(m/n)

Author(s)

Roger Koenker

References

Bose, A. and S. Chatterjee, (2003) Generalized bootstrap for estimators of minimizers of convex functions, J. Stat. Planning and Inf, 117, 225-239. Portnoy, S. (2013) The Jackknife's Edge: Inference for Censored Quantile Regression, CSDA, forthcoming.

See Also

summary.crq

Description

These functions can be used to construct standard errors, confidence intervals and tests of hypotheses regarding quantile regression models.

Usage

boot.rq(x, y, tau = 0.5, R = 200, bsmethod = "xy", mofn = length(y), 
	coef = NULL, blbn = NULL, cluster = NULL, U = NULL,  ...)

Arguments

Details

Their are several refinements that are still unimplemented. Percentile methods should be incorporated, and extensions of the methods to be used in anova.rq should be made. And more flexibility about what algorithm is used would also be good.

Value

A list consisting of two elements: A matrix B of dimension R by p is returned with the R resampled estimates of the vector of quantile regression parameters. When mofn < n for the "xy" method this matrix has been deflated by the factor sqrt(m/n). A matrix U of sampled indices (for bsmethod in c("xy", "wxy")) or gradient evaluations (for bsmethod in c("pwy", "cluster")) used to generate the bootstrapped realization, and potentially reused for other taus when invoked from summary.rqs.

Author(s)

Roger Koenker (and Xuming He and M. Kocherginsky for the mcmb code)

References

[1] Koenker, R. W. (1994). Confidence Intervals for regression quantiles, in P. Mandl and M. Huskova (eds.), Asymptotic Statistics, 349–359, Springer-Verlag, New York.

[2] Kocherginsky, M., He, X. and Mu, Y. (2005). Practical Confidence Intervals for Regression Quantiles, Journal of Computational and Graphical Statistics, 14, 41-55.

[3] Hagemann, A. (2017) Cluster Robust Bootstrap inference in quantile regression models, Journal of the American Statistical Association , 112, 446–456.

[4] He, X. and Hu, F. (2002). Markov Chain Marginal Bootstrap. Journal of the American Statistical Association , Vol. 97, no. 459, 783-795.

[5] Parzen, M. I., L. Wei, and Z. Ying (1994): A resampling method based on pivotal estimating functions,” Biometrika, 81, 341–350.

[6] Bose, A. and S. Chatterjee, (2003) Generalized bootstrap for estimators of minimizers of convex functions, J. Stat. Planning and Inf, 117, 225-239.

[7] Chamberlain G. and Imbens G.W. (2003) Nonparametric Applications of Bayesian Inference, Journal of Business & Economic Statistics, 21, pp. 12-18.

[8] Feng, Xingdong, Xuming He, and Jianhua Hu (2011) Wild Bootstrap for Quantile Regression, Biometrika, 98, 995–999.

See Also

summary.rq

Examples

y <- rnorm(50)
x <- matrix(rnorm(100),50)
fit <- rq(y~x,tau = .4)
summary(fit,se = "boot", bsmethod= "xy")
summary(fit,se = "boot", bsmethod= "pwy")
#summary(fit,se = "boot", bsmethod= "mcmb")

Description

Bootstrap method exploiting preprocessing strategy to reduce computation time for large problem. In contrast to boot.rq.pxy which uses the classical multinomial sampling scheme and is coded in R, this uses the exponentially weighted bootstrap scheme and is coded in fortran and consequently is considerably faster in larger problems.

Usage

boot.rq.pwxy(x, y, tau, coef, R = 200, m0 = NULL, eps = 1e-06, ...)

Arguments

Details

The fortran implementation is quite similar to the R code for boot.rq.pxy except that there is no multinomial sampling. Instead rexp(n) weights are used.

Value

returns a list with elements:

Author(s)

Blaise Melly and Roger Koenker

References

Chernozhukov, V. I. Fernandez-Val and B. Melly, Fast Algorithms for the Quantile Regression Process, 2019, arXiv, 1909.05782,

Portnoy, S. and R. Koenker, The Gaussian Hare and the Laplacian Tortoise, Statistical Science, (1997) 279-300

See Also

boot.rq.pxy

Description

Bootstrap method exploiting preprocessing strategy to reduce computation time for large problem.

Usage

boot.rq.pxy(x, y, s, tau = 0.5, coef, method = "fn", Mm.factor = 3)

Arguments

Details

See references for further details.

Value

Returns matrix of bootstrap estimates.

Author(s)

Blaise Melly and Roger Koenker

References

Chernozhukov, V. I. Fernandez-Val and B. Melly, Fast Algorithms for the Quantile Regression Process, 2019, arXiv, 1909.05782,

Portnoy, S. and R. Koenker, The Gaussian Hare and the Laplacian Tortoise, Statistical Science, (1997) 279-300

See Also

rq.fit.ppro

Description

Boscovich data used to estimate the ellipticity of the earth. There are five measurements of the arc length of one degree of latitude taken at 5 different latitudes. See Koenker (2005) for further details and references.

Usage

data(Bosco)

Format

A data frame containing 5 observations on 2 variables

x

sine squared of latitude measured in degrees

y

arc length of one degree of latitude measured in toise - 56,700, one toise approximately equals 1.95 meters.

References

Koenker, R. (2005), "Quantile Regression", Cambridge.

Examples

data(Bosco)
plot(0:10/10,0:10*100,xlab="sin^2(latitude)",
        ylab="arc-length of 1 degree of latitude",type="n")
points(Bosco)
text(Bosco, pos = 3, rownames(Bosco))
z <- rq(y ~ x, tau = -1, data = Bosco)
title("Boscovitch Ellipticity of the Earth Example")
xb <- c(.85,.9,.6,.6)
yb <- c(400,600,450,600)
for(i in 1:4){
        abline(c(z$sol[4:5,i]))
        interval <- paste("t=(",format(round(z$sol[1,i],2)),",",
                format(round(z$sol[1,i+1],2)),")",delim="")
        text(xb[i],yb[i],interval)
        }

Description

Cobar Ore data from Green and Silverman (1994). The data consists of measurements on the "true width" of an ore-bearing rock layer from a mine in Cobar, Australia.

Usage

data(CobarOre)

Format

A data frame with 38 observations on the following 3 variables.

x

x-coordinate of location of mine site

y

y-coordinate of location of mine site

z

ore thickness

Source

Green, P.J. and B.W. Silverman (1994) Nonparametric Regression Generalized Linear Models: A roughness penalty approach, Chapman Hall.

Examples

data(CobarOre)
plot(CobarOre)

Description

All m combinations of the first n integers taken p at a time are computed and return as an p by m matrix. The columns of the matrix are ordered so that adjacent columns differ by only one element. This is just a reordered version of combn in base R, but the ordering is useful for some applications.

Usage

combos(n,p)

Arguments

Value

a matrix of dimension p by choose(n,p)

Note

Implementation based on a Pascal algorithm of Limin Xiang and Kazuo Ushijima (2001) translated to ratfor for R. If you have rgl installed you might try demo("combos") for a visual impression of how this works.

References

Limin Xiang and Kazuo Ushijima (2001) "On O(1) Time Algorithms for Combinatorial Generation," Computer Journal, 44(4), 292-302.

Examples

H <- combos(20,3)

Description

Critical values for uniform confidence bands for rqss fitting

Usage

critval(kappa, alpha = 0.05, rdf = 0)

Arguments

Details

The Hotelling tube approach to inference has a long and illustrious history. See Johansen and Johnstone (1989) for an overview. The implementation here is based on Sun and Loader (1994) and Loader's locfit package, although a simpler root finding approach is substituted for the iterative method used there. At this stage, only univariate bands may be constructed.

Value

A scalar critical value that acts as a multiplier for the uniform confidence band construction.

References

Hotelling, H. (1939): “Tubes and Spheres in $n$-spaces, and a class of statistical problems,” Am J. Math, 61, 440–460.

Johansen, S., I.M. Johnstone (1990): “Hotelling's Theorem on the Volume of Tubes: Some Illustrations in Simultaneous Inference and Data Analysis,” The Annals of Statistics, 18, 652–684.

Sun, J. and C.V. Loader: (1994) “Simultaneous Confidence Bands for Linear Regression and smoothing,” The Annals of Statistics, 22, 1328–1345.

See Also

plot.rqss

Description

Fits a conditional quantile regression model for censored data. There are three distinct methods: the first is the fixed censoring method of Powell (1986) as implemented by Fitzenberger (1996), the second is the random censoring method of Portnoy (2003). The third method is based on Peng and Huang (2008).

Usage

crq(formula, taus, data, subset, weights, na.action, 
	method = c("Powell", "Portnoy", "Portnoy2", "PengHuang"), contrasts = NULL, ...)
crq.fit.pow(x, y, yc, tau=0.5, weights=NULL, start, left=TRUE, maxit = 500)
crq.fit.pen(x, y, cen, weights=NULL, grid, ctype = "right") 
crq.fit.por(x, y, cen, weights=NULL, grid, ctype = "right") 
crq.fit.por2(x, y, cen, weights=NULL, grid, ctype = "right") 
Curv(y, yc, ctype=c("left","right"))
## S3 method for class 'crq'
print(x, ...)
## S3 method for class 'crq'
print(x, ...)
## S3 method for class 'crq'
predict(object, newdata,  ...)
## S3 method for class 'crqs'
predict(object, newdata, type = NULL, ...)
## S3 method for class 'crq'
coef(object,taus = 1:4/5,...)

Arguments

Details

The Fitzenberger algorithm uses a variant of the Barrodale and Roberts simplex method. Exploiting the fact that the solution must be characterized by an exact fit to p points when there are p parameters to be estimated, at any trial basic solution it computes the directional derivatives in the 2p distinct directions and choses the direction that (locally) gives steepest descent. It then performs a one-dimensional line search to choose the new basic observation and continues until it reaches a local mimumum. By default it starts at the naive rq solution ignoring the censoring; this has the (slight) advantage that the estimator is consequently equivariant to canonical transformations of the data. Since the objective function is no longer convex there can be no guarantee that this produces a global minimum estimate. In small problems exhaustive search over solutions defined by p-element subsets of the n observations can be used, but this quickly becomes impractical for large p and n. This global version of the Powell estimator can be invoked by specifying start = "global". Users interested in this option would be well advised to compute choose(n,p) for their problems before trying it. The method operates by pivoting through this many distinct solutions and choosing the one that gives the minimal Powell objective. The algorithm used for the Portnoy method is described in considerable detail in Portnoy (2003). There is a somewhat simplified version of the Portnoy method that is written in R and iterates over a discrete grid. This version should be considered somewhat experimental at this stage, but it is known to avoid some difficulties with the more complicated fortran version of the algorithm that can occur in degenerate problems. Both the Portnoy and Peng-Huang estimators may be unable to compute estimates of the conditional quantile parameters in the upper tail of distribution. Like the Kaplan-Meier estimator, when censoring is heavy in the upper tail the estimated distribution is defective and quantiles are only estimable on a sub-interval of (0,1). The Peng and Huang estimator can be viewed as a generalization of the Nelson Aalen estimator of the cumulative hazard function, and can be formulated as a variant of the conventional quantile regression dual problem. See Koenker (2008) for further details. This paper is available from the package with vignette("crq").

Value

An object of class crq.

Author(s)

Steve Portnoy and Roger Koenker

References

Fitzenberger, B. (1996): “A Guide to Censored Quantile Regressions,” in Handbook of Statistics, ed. by C.~Rao, and G.~Maddala. North-Holland: New York.

Fitzenberger, B. and P. Winker (2007): “Improving the Computation of Censored Quantile Regression Estimators,” CSDA, 52, 88-108.

Koenker, R. (2008): “Censored Quantile Regression Redux,” J. Statistical Software, 27, https://www.jstatsoft.org/v27/i06.

Peng, L and Y Huang, (2008) Survival Analysis with Quantile Regression Models, J. Am. Stat. Assoc., 103, 637-649.

Portnoy, S. (2003) “Censored Quantile Regression,” JASA, 98,1001-1012.

Powell, J. (1986) “Censored Regression Quantiles,” J. Econometrics, 32, 143–155.

See Also

summary.crq

Examples

# An artificial Powell example
set.seed(2345)
x <- sqrt(rnorm(100)^2)
y <-  -0.5 + x +(.25 + .25*x)*rnorm(100)
plot(x,y, type="n")
s <- (y > 0)
points(x[s],y[s],cex=.9,pch=16)
points(x[!s],y[!s],cex=.9,pch=1)
yLatent <- y
y <- pmax(0,y)
yc <- rep(0,100)
for(tau in (1:4)/5){
        f <- crq(Curv(y,yc) ~ x, tau = tau, method = "Pow")
        xs <- sort(x)
        lines(xs,pmax(0,cbind(1,xs)%*%f$coef),col="red")
        abline(rq(y ~ x, tau = tau), col="blue")
        abline(rq(yLatent ~ x, tau = tau), col="green")
        }
legend(.15,2.5,c("Naive QR","Censored QR","Omniscient QR"),
        lty=rep(1,3),col=c("blue","red","green"))

# crq example with left censoring
set.seed(1968)
n <- 200
x <-rnorm(n)
y <- 5 + x + rnorm(n)
plot(x,y,cex = .5)
c <- 4 + x + rnorm(n)
d <- (y > c)
points(x[!d],y[!d],cex = .5, col = 2)
f <- crq(survival::Surv(pmax(y,c), d, type = "left") ~ x, method = "Portnoy")
g <- summary(f)
for(i in 1:4) abline(coef(g[[i]])[,1])

Description

With malice aforethought, dither adds a specified random perturbation to each element of the input vector, usually employed as a device to mitigate the effect of ties.

Usage

dither(x, type = "symmetric", value = NULL)

Arguments

Details

The function dither operates slightly differently than the function jitter in base R, permitting strictly positive perturbations with the option type = "right" and using somewhat different default schemes for the scale of the perturbation. Dithering the response variable is frequently a useful option in quantile regression fitting to avoid deleterious effects of degenerate solutions. See, e.g. Machado and Santos Silva (2005). For a general introduction and some etymology see the Wikipedia article on "dither". For integer data it is usually advisable to use value = 1. When 'x' is a matrix or array dither treats all elements as a vector but returns an object of the original class.

Value

A dithered version of the input vector 'x'.

Note

Some further generality might be nice, for example something other than uniform noise would be desirable in some circumstances. Note that when dithering you are entering into the "state of sin" that John von Neumann famously attributed to anyone considering "arithmetical methods of producing random digits." If you need to preserve reproducibility, then set.seed is your friend.

Author(s)

R. Koenker

References

Machado, J.A.F. and Santos Silva, J.M.C. (2005), Quantiles for Counts, Journal of the American Statistical Association, vol. 100, no. 472, pp. 1226-1237.

See Also

jitter

Examples

x <- rlnorm(40)
y <- rpois(40, exp(.5 + log(x)))
f <- rq(dither(y, type = "right", value = 1) ~ x)

Description

Interface to rq.fit and rq.wfit for fitting dynamic linear quantile regression models. The interface is based very closely on Achim Zeileis's dynlm package. In effect, this is mainly “syntactic sugar” for formula processing, but one should never underestimate the value of good, natural sweeteners.

Usage

dynrq(formula, tau = 0.5, data, subset, weights, na.action, method = "br",
  contrasts = NULL, start = NULL, end = NULL, ...)

Arguments

Details

The interface and internals of dynrq are very similar to rq, but currently dynrq offers two advantages over the direct use of rq for time series applications of quantile regression: extended formula processing, and preservation of time series attributes. Both features have been shamelessly lifted from Achim Zeileis's package dynlm.

For specifying the formula of the model to be fitted, there are several functions available which allow for convenient specification of dynamics (via d() and L()) or linear/cyclical patterns (via trend(), season(), and harmon()). These new formula functions require that their arguments are time series objects (i.e., "ts" or "zoo").

Dynamic models: An example would be d(y) ~ L(y, 2), where d(x, k) is diff(x, lag = k) and L(x, k) is lag(x, lag = -k), note the difference in sign. The default for k is in both cases 1. For L(), it can also be vector-valued, e.g., y ~ L(y, 1:4).

Trends: y ~ trend(y) specifies a linear time trend where (1:n)/freq is used by default as the covariate, n is the number of observations and freq is the frequency of the series (if any, otherwise freq = 1). Alternatively, trend(y, scale = FALSE) would employ 1:n and time(y) would employ the original time index.

Seasonal/cyclical patterns: Seasonal patterns can be specified via season(x, ref = NULL) and harmonic patterns via harmon(x, order = 1). season(x, ref = NULL) creates a factor with levels for each cycle of the season. Using the ref argument, the reference level can be changed from the default first level to any other. harmon(x, order = 1) creates a matrix of regressors corresponding to cos(2 * o * pi * time(x)) and sin(2 * o * pi * time(x)) where o is chosen from 1:order.

See below for examples.

Another aim of dynrq is to preserve time series properties of the data. Explicit support is currently available for "ts" and "zoo" series. Internally, the data is kept as a "zoo" series and coerced back to "ts" if the original dependent variable was of that class (and no internal NAs were created by the na.action).

See Also

zoo, merge.zoo

Examples

###########################
## Dynamic Linear Quantile Regression Models ##
###########################

if(require(zoo)){
## multiplicative median SARIMA(1,0,0)(1,0,0)_12 model fitted to UK seatbelt data
     uk <- log10(UKDriverDeaths)
     dfm <- dynrq(uk ~ L(uk, 1) + L(uk, 12))
     dfm

     dfm3 <- dynrq(uk ~ L(uk, 1) + L(uk, 12),tau = 1:3/4)
     summary(dfm3)
 ## explicitly set start and end
     dfm1 <- dynrq(uk ~ L(uk, 1) + L(uk, 12), start = c(1975, 1), end = c(1982, 12))
 ## remove lag 12
     dfm0 <- update(dfm1, . ~ . - L(uk, 12))
     tuk1  <- anova(dfm0, dfm1)
 ## add seasonal term
     dfm1 <- dynrq(uk ~ 1, start = c(1975, 1), end = c(1982, 12))
     dfm2 <- dynrq(uk ~ season(uk), start = c(1975, 1), end = c(1982, 12))
     tuk2 <- anova(dfm1, dfm2)
 ## regression on multiple lags in a single L() call
     dfm3 <- dynrq(uk ~ L(uk, c(1, 11, 12)), start = c(1975, 1), end = c(1982, 12))
     anova(dfm1, dfm3)
}

###############################
## Time Series Decomposition ##
###############################

## airline data
## Not run: 
ap <- log(AirPassengers)
fm <- dynrq(ap ~ trend(ap) + season(ap), tau = 1:4/5)
sfm <- summary(fm)
plot(sfm)

## End(Not run)

## Alternative time trend specifications:
##   time(ap)                  1949 + (0, 1, ..., 143)/12
##   trend(ap)                 (1, 2, ..., 144)/12
##   trend(ap, scale = FALSE)  (1, 2, ..., 144)

###############################
## An Edgeworth (1886) Problem##
###############################
# DGP
## Not run: 
fye <- function(n, m = 20){
    a <- rep(0,n)
    s <- sample(0:9, m, replace = TRUE)
    a[1] <- sum(s)
    for(i in 2:n){
       s[sample(1:20,1)] <- sample(0:9,1)
       a[i] <- sum(s)
    }
    zoo::zoo(a)
}
x <- fye(1000)
f <- dynrq(x ~ L(x,1))
plot(x,cex = .5, col = "red")
lines(fitted(f), col = "blue")

## End(Not run)

Description

Engel food expenditure data used in Koenker and Bassett(1982). This is a regression data set consisting of 235 observations on income and expenditure on food for Belgian working class households.

Usage

data(engel)

Format

A data frame containing 235 observations on 2 variables

income

annual household income in Belgian francs

foodexp

annual household food expenditure in Belgian francs

References

Koenker, R. and Bassett, G (1982) Robust Tests of Heteroscedasticity based on Regression Quantiles; Econometrica 50, 43–61.

Examples

## See also    demo("engel1")
##             --------------

data(engel)
plot(engel, log = "xy",
     main = "'engel' data  (log - log scale)")
plot(log10(foodexp) ~ log10(income), data = engel,
     main = "'engel' data  (log10 - transformed)")
taus <- c(.15, .25, .50, .75, .95, .99)
rqs <- as.list(taus)
for(i in seq(along = taus)) {
  rqs[[i]] <- rq(log10(foodexp) ~ log10(income), tau = taus[i], data = engel)
  lines(log10(engel$income), fitted(rqs[[i]]), col = i+1)
}
legend("bottomright", paste("tau = ", taus), inset = .04,
       col = 2:(length(taus)+1), lty=1)

Description

Show the FAQ or ChangeLog of a specified package

Usage

FAQ(pkg = "quantreg")
ChangeLog(pkg = "quantreg")

Arguments

Details

Assumes that the FAQ and/or ChangeLog files exist in the proper "inst" directory.

Value

Has only the side effect of showing the files on the screen.

Description

Time Series of Weekly US Gasoline Prices: 1990:8 – 2003:26

Usage

data("gasprice")

Examples

data(gasprice)

Description

Tests of the hypothesis that a linear model specification is of the location shift or location-scale shift form. The tests are based on the Doob-Meyer Martingale transformation approach proposed by Khmaladze(1981) for general goodness of fit problems as adapted to quantile regression by Koenker and Xiao (2002).

Usage

KhmaladzeTest(formula, data = NULL, taus = 1:99/100, nullH = "location" ,  
	trim = c(0.05, 0.95), h = 1, ...)

Arguments

Value

an object of class KhmaladzeTest is returned containing:

References

Khmaladze, E. (1981) “Martingale Approach in the Theory of Goodness-of-fit Tests,” Theory of Prob. and its Apps, 26, 240–257.

Koenker, Roger and Zhijie Xiao (2002), “Inference on the Quantile Regression Process”, Econometrica, 81, 1583–1612. http://www.econ.uiuc.edu/~roger/research/inference/inference.html

Examples

data(barro)
T = KhmaladzeTest( y.net ~ lgdp2 + fse2 + gedy2 + Iy2 + gcony2, 
		data = barro, taus = seq(.05,.95,by = .01))
plot(T)

Description

The function 'kuantile' computes sample quantiles corresponding to the specified probabilities. The intent is to mimic the generic (base) function 'quantile' but using a variant of the Floyd and Rivest (1975) algorithm which is somewhat quicker, especially for large sample sizes.

Usage

kuantile(x, probs = seq(0, 1, .25), na.rm = FALSE, names = TRUE, type = 7, ...)

Arguments

Details

A vector of length 'length(p)' is returned. See the documentation for 'quantile' for further details on the types. The algorithm was written by K.C. Kiwiel. It is a modified version of the (algol 68) SELECT procedure of Floyd and Rivest (1975), incorporating modifications of Brown(1976). The algorithm has linear growth in the number of comparisons required as sample size grows. For the median, average case behavior requires $1.5 n + O((n log n)^{1/2})$ comparisons. See Kiwiel (2005) and Knuth (1998) for further details. When the number of required elements of p is large, it may be preferable to revert to a full sort.

Value

A vector of quantiles of the same length as the vector p.

Author(s)

K.C. Kiwiel, R interface: Roger Koenker

References

R.W. Floyd and R.L. Rivest: "Algorithm 489: The Algorithm SELECT—for Finding the $i$th Smallest of $n$ Elements", Comm. ACM 18, 3 (1975) 173,

T. Brown: "Remark on Algorithm 489", ACM Trans. Math. Software 3, 2 (1976), 301-304.

K.C. Kiwiel: On Floyd and Rivest's SELECT Algorithm, Theoretical Computer Sci. 347 (2005) 214-238.

D. Knuth, The Art of Computer Programming, Volume 3, Sorting and Searching, 2nd Ed., (1998), Addison-Wesley.

See Also

quantile

Examples

kuantile(x <- rnorm(1001))# Extremes & Quartiles by default

     ### Compare different types
     p <- c(0.1,0.5,1,2,5,10,50)/100
     res <- matrix(as.numeric(NA), 9, 7)
     for(type in 1:9) res[type, ] <- y <- kuantile(x,  p, type=type)
     dimnames(res) <- list(1:9, names(y))
     ktiles <- res

     ### Compare different types
     p <- c(0.1,0.5,1,2,5,10,50)/100
     res <- matrix(as.numeric(NA), 9, 7)
     for(type in 1:9) res[type, ] <- y <- quantile(x,  p, type=type)
     dimnames(res) <- list(1:9, names(y))
     qtiles <- res

     max(abs(ktiles - qtiles))

Description

Default procedure for selection of lambda in lasso constrained quantile regression as proposed by Belloni and Chernozhukov (2011)

Usage

LassoLambdaHat(X, R = 1000, tau = 0.5, C = 1, alpha = 0.95)

Arguments

Details

As proposed by Belloni and Chernozhukov, a reasonable default lambda would be the upper quantile of the simulated values. The procedure is based on idea that a simulated gradient can be used as a pivotal statistic. Elements of the default vector are standardized by the respective standard deviations of the covariates. Note that the sqrt(tau(1-tau)) factor cancels in their (2.4) (2.6). In this formulation even the intercept is penalized. If the lower limit of the simulated interval is desired one can specify alpha = 0.05.

Value

vector of default lambda values of length p, the column dimension of X.

References

Belloni, A. and V. Chernozhukov. (2011) l1-penalized quantile regression in high-dimensional sparse models. Annals of Statistics, 39 82 - 130.

Examples

n <- 200
p <- 10
x <- matrix(rnorm(n*p), n, p)
b <- c(1,1, rep(0, p-2))
y <- x %*% b + rnorm(n)
f <- rq(y ~ x, tau = 0.8, method = "lasso")
# See f$lambda to see the default lambda selection

Description

Generic function for converting an R object into a latex file.

Usage

latex(x, ...)

Arguments

See Also

latex.table, latex.summary.rqs

Description

Produces a file with latex commands for a table of rq results.

Usage

## S3 method for class 'summary.rqs'
latex(x, transpose = FALSE, caption = "caption goes here.", 
		digits = 3, file = as.character(substitute(x)), ...)

Arguments

Details

Calls latex.table.

Value

Returns invisibly after writing the file.

Author(s)

R. Koenker

See Also

summary.rqs, latex.table

Description

Automatically generates a latex formatted table from the matrix x Controls rounding, alignment, etc, etc

Usage

## S3 method for class 'table'
latex(x, file=as.character(substitute(x)), 
	rowlabel=file, rowlabel.just="l", cgroup, n.cgroup, rgroup, n.rgroup=NULL, 
	digits, dec, rdec, cdec, append=FALSE, dcolumn=FALSE, cdot=FALSE, 
	longtable=FALSE, table.env=TRUE, lines.page=40, caption, caption.lot, 
	label=file, double.slash=FALSE,...)

Arguments

Value

returns invisibly

Author(s)

Roger Koenker

References

Minor modification of Frank Harrell's Splus code

Description

This function fits a linear model by recursive least squares. It is a utility routine for the KhmaladzeTest function of the quantile regression package.

Usage

lm.fit.recursive(X, y, int=TRUE)

Arguments

Value

return p by n matrix of fitted parameters, where p. The ith column gives the solution up to "time" i.

Author(s)

R. Koenker

References

A. Harvey, (1993) Time Series Models, MIT

Description

This is a toy function to illustrate how to do locally polynomial quantile regression univariate smoothing.

Usage

lprq(x, y, h, tau = .5, m = 50)

Arguments

Details

The function obviously only does locally linear fitting but can be easily adapted to locally polynomial fitting of higher order. The author doesn't really approve of this sort of smoothing, being more of a spline person, so the code is left is its (almost) most trivial form.

Value

The function compute a locally weighted linear quantile regression fit at each of the m design points, and returns:

Note

One can also consider using B-spline expansions see bs.

Author(s)

R. Koenker

References

Koenker, R. (2004) Quantile Regression

See Also

rqss for a general approach to oonparametric QR fitting.

Examples

require(MASS)
data(mcycle)
attach(mcycle)
plot(times,accel,xlab = "milliseconds", ylab = "acceleration (in g)")
hs <- c(1,2,3,4)
for(i in hs){
        h = hs[i]
        fit <- lprq(times,accel,h=h,tau=.5)
        lines(fit$xx,fit$fv,lty=i)
        }
legend(50,-70,c("h=1","h=2","h=3","h=4"),lty=1:length(hs))

Description

Observations on the maximal running speed of mammal species and their body mass.

Usage

data(Mammals)

Format

A data frame with 107 observations on the following 4 variables.

weight

Body mass in Kg for "typical adult sizes"

speed

Maximal running speed (fastest sprint velocity on record)

hoppers

logical variable indicating animals that ambulate by hopping, e.g. kangaroos

specials

logical variable indicating special animals with "lifestyles in which speed does not figure as an important factor": Hippopotamus, raccoon (Procyon), badger (Meles), coati (Nasua), skunk (Mephitis), man (Homo), porcupine (Erithizon), oppossum (didelphis), and sloth (Bradypus)

Details

Used by Chappell (1989) and Koenker, Ng and Portnoy (1994) to illustrate the fitting of piecewise linear curves.

Source

Garland, T. (1983) The relation between maximal running speed and body mass in terrestrial mammals, J. Zoology, 199, 1557-1570.

References

Koenker, R., P. Ng and S. Portnoy, (1994) Quantile Smoothing Splines” Biometrika, 81, 673-680.

Chappell, R. (1989) Fitting Bent Lines to Data, with Applications ot Allometry, J. Theo. Biology, 138, 235-256.

See Also

rqss

Examples

data(Mammals)
attach(Mammals)
x <- log(weight)
y <- log(speed)
plot(x,y, xlab="Weight in log(Kg)", ylab="Speed in log(Km/hour)",type="n")
points(x[hoppers],y[hoppers],pch = "h", col="red")
points(x[specials],y[specials],pch = "s", col="blue")
others <- (!hoppers & !specials)
points(x[others],y[others], col="black",cex = .75)
fit <- rqss(y ~ qss(x, lambda = 1),tau = .9)
plot(fit)

Description

Daily maximum temperatures in Melbourne, Australia, from 1981-1990. Leap days have been omitted.

Usage

data(MelTemp)

Format

Time series of frequency 365

Source

Hyndman, R.J., Bashtannyk, D.M. and Grunwald, G.K. (1996) "Estimating and visualizing conditional densities". _Journal of Computational and Graphical Statistics_, *5*, 315-336.

Examples

data(MelTemp)
demo(Mel)

Description

function to recursively substitute arguments into rqss formula

Usage

Munge(formula, ...)

Arguments

Details

Intended (originally) for use with demo(MCV). Based on an R-help suggestion of Gabor Grothendieck.

Value

A new formula after substitution

See Also

demo(MCV)

Examples

lams <- c(1.3, 3.3)
f <- y ~ qss(x, lambda = lams[1]) + qss(z, lambda = lams[2]) + s
ff <- Munge(f, lams = lams)

Description

This function implements an R version of an interior point method for computing the solution to quantile regression problems which are nonlinear in the parameters. The algorithm is based on interior point ideas described in Koenker and Park (1994).

Usage

nlrq(formula, data=parent.frame(), start, tau=0.5, 
	control, trace=FALSE,method="L-BFGS-B")
## S3 method for class 'nlrq'
summary(object, ...)
## S3 method for class 'summary.nlrq'
print(x, digits = max(5, .Options$digits - 2), ...)

Arguments

Details

An ‘nlrq’ object is a type of fitted model object. It has methods for the generic functions ‘coef’ (parameters estimation at best solution), ‘formula’ (model used), ‘deviance’ (value of the objective function at best solution), ‘print’, ‘summary’, ‘fitted’ (vector of fitted variable according to the model), ‘predict’ (vector of data points predicted by the model, using a different matrix for the independent variables) and also for the function ‘tau’ (quantile used for fitting the model, as the tau argument of the function). Further help is also available for the method ‘residuals’. The summary method for nlrq uses a bootstrap approach based on the final linearization of the model evaluated at the estimated parameters.

Value

A list consisting of:

Author(s)

Based on S code by Roger Koenker modified for R and to accept models as specified by nls by Philippe Grosjean.

References

Koenker, R. and Park, B.J. (1994). An Interior Point Algorithm for Nonlinear Quantile Regression, Journal of Econometrics, 71(1-2): 265-283.

See Also

nlrq.control , residuals.nlrq

Examples

# build artificial data with multiplicative error
Dat <- NULL; Dat$x <- rep(1:25, 20)
set.seed(1)
Dat$y <- SSlogis(Dat$x, 10, 12, 2)*rnorm(500, 1, 0.1)
plot(Dat)
# fit first a nonlinear least-square regression
Dat.nls <- nls(y ~ SSlogis(x, Asym, mid, scal), data=Dat); Dat.nls
lines(1:25, predict(Dat.nls, newdata=list(x=1:25)), col=1)
# then fit the median using nlrq
Dat.nlrq <- nlrq(y ~ SSlogis(x, Asym, mid, scal), data=Dat, tau=0.5, trace=TRUE)
lines(1:25, predict(Dat.nlrq, newdata=list(x=1:25)), col=2)
# the 1st and 3rd quartiles regressions
Dat.nlrq <- nlrq(y ~ SSlogis(x, Asym, mid, scal), data=Dat, tau=0.25, trace=TRUE)
lines(1:25, predict(Dat.nlrq, newdata=list(x=1:25)), col=3)
Dat.nlrq <- nlrq(y ~ SSlogis(x, Asym, mid, scal), data=Dat, tau=0.75, trace=TRUE)
lines(1:25, predict(Dat.nlrq, newdata=list(x=1:25)), col=3)
# and finally "external envelopes" holding 95 percent of the data
Dat.nlrq <- nlrq(y ~ SSlogis(x, Asym, mid, scal), data=Dat, tau=0.025, trace=TRUE)
lines(1:25, predict(Dat.nlrq, newdata=list(x=1:25)), col=4)
Dat.nlrq <- nlrq(y ~ SSlogis(x, Asym, mid, scal), data=Dat, tau=0.975, trace=TRUE)
lines(1:25, predict(Dat.nlrq, newdata=list(x=1:25)), col=4)
leg <- c("least squares","median (0.5)","quartiles (0.25/0.75)",".95 band (0.025/0.975)")
legend(1, 12.5, legend=leg, lty=1, col=1:4)

Description

Set algorithmic parameters for nlrq (nonlinear quantile regression function)

Usage

nlrq.control(maxiter=100, k=2, InitialStepSize = 1, big=1e+20, eps=1e-07, beta=0.97)

Arguments

See Also

nlrq

Description

Estimation and inference about the tail behavior of the response in linear models are based on the adaptation of the univariate Hill (1975) and Pickands (1975) estimators for quantile regression by Chernozhukov, Fernandez-Val and Kaji (2018).

Usage

ParetoTest(formula, tau = 0.1, data = NULL, flavor = "Hill", m = 2, cicov = .9, ...)

Arguments

Value

an object of class ParetoTest is returned containing:

References

Chernozhukov, Victor, Ivan Fernandez-Val, and Tetsuya Kaji, (2018) Extremal Quantile Regression, in Handbook of Quantile Regression, Eds. Roger Koenker, Victor Chernozhukov, Xuming He, Limin Peng, CRC Press.

Hill, B. M. (1975). A simple general approach to inference about the tail of a distribution. The Annals of Statistics 3(5), 1163-1174.

Pickands, J. (1975). Statistical inference using extreme order statistics. The Annals of Statistics 3(1), 119-131.

Examples

n = 500
x = rnorm(n)
y = x + rt(n,2)
Z = ParetoTest(y ~ x, .9, flavor = "Pickands")

Description

Data from sequence experiments conducted by C.S. Pierce in 1872 to determine the distribution of response times to an auditory stimulus.

Usage

data(Peirce)

Format

A link{list} of 24 objects each representing one day of the experiment. Each element of the list consists of three components: the date the measurements were made, an x component recording the response time in milliseconds, and an associated y component recording a count of the number of times that the response was recorded to be equal to be equal to the corresponding x entry. There are roughly 500 observations (counts) on each of the 24 days.

Details

A detailed description of the experiment can be found in Peirce (1873). A young man of about 18 with no prior experience was employed to respond to a signal “consisting of a sharp sound like a rap, the answer being made upon a telegraph-operator's key nicely adjusted.” The response times, made with the aid of a Hipp cronoscope were recorded to the nearest millisecond. The data was analyzed by Peirce who concluded that after the first day, when the the observer was entirely inexperienced, the curves representing the densities of the response times “differed very little from that derived from the theory of least squares,” i.e. from the Gaussian density.

The data was subsequently analysed by Samama, in a diploma thesis supervised by Maurice Frechet, who reported briefly the findings in Frechet (1924), and by Wilson and Hilferty (1929). In both instances the reanalysis showed that Laplace's first law of error, the double exponential distribution, was a better representation for the data than was the Gaussian law. Koenker (2009) constains further discussion and an attempt to reproduce the Wilson and Hilferty analysis.

The data is available in two formats: The first in a "raw" form as 24 text files as scanned from the reprinted Peirce source, the second as an R dataset Peirce.rda containing the list. Only the latter is provided here, for the raw data and how to read see the more complete archive at: http://www.econ.uiuc.edu/~roger/research/frechet/frechet.html See the examples section below for some details on provisional attempt to reproduce part of the Wilson and Hilferty analysis. An open question regarding the dataset is: How did Wilson and Hilferty compute standard deviations for the median as they appear in their table? The standard textbook suggestion of Yule (1917) yields far too small a bandwidth. The methods employed in the example section below, which rely on relatively recent proposals, are somewhat closer, but still deviate somewhat from the results reported by Wilson and Hilferty.

Source

Peirce, C.~S. (1873): “On the Theory of Errors of Observation,” Report of the Superintendent of the U.S. Coast Survey, pp. 200–224, Reprinted in The New Elements of Mathematics, (1976) collected papers of C.S. Peirce, ed. by C. Eisele, Humanities Press: Atlantic Highlands, N.J., vol. 3, part 1, 639–676.

References

Fr\'echet, M. (1924): “Sur la loi des erreurs d'observation,” Matematichiskii Sbornik, 32, 5–8. Koenker, R. (2009): “The Median is the Message: Wilson and Hilferty's Reanalysis of C.S. Peirce's Experiments on the Law of Errors,” American Statistician, 63, 20-25. Wilson, E.~B., and M.~M. Hilferty (1929): “Note on C.S. Peirces Experimental Discussion of the Law of Errors,” Proceedings of the National Academy of Sciences of the U.S.A., 15, 120–125. Yule, G.~U. (1917): An Introduction to the Theory of Statistics. Charles Griffen: London, 4 edn.

Examples

# Make table like Wilson and Hilferty

data("Peirce")
set.seed(10) #Dither the counts
tab <- matrix(0,24,11)
for(i in 1:24){
	y <- rep(Peirce[[i]]$x, Peirce[[i]]$y) + runif(sum(Peirce[[i]]$y), -.5, .5)
	f1 <- summary(rq(y~1),se="iid")$coef[1:2]
	n <- length(y)
	f0 <- 1/(2 * sum(abs(y-f1[1])/n)) #Laplace proposal
	f0 <- (1/(2 * f0))/ sqrt(n)
	f2 <- summary(lm(y~1))$coef[1:2]
	outm <- sum(y < (f1[1] - 3.1 * sqrt(n) * f2[2]))
	outp <- sum(y > (f1[1] + 3.1 * sqrt(n) * f2[2]))
	outt <- outm + outp
	inm <- y > (f1[1] - 0.25 * sqrt(n) * f2[2])
	inp <- y < (f1[1] + 0.25 * sqrt(n) * f2[2])
	int <- sum(inm * inp)
	Eint <- round(n * (pnorm(.25) - pnorm(-.25)))
	excess <- round(100*(int - Eint)/Eint)
	tab[i,] <- c(f1, f0, f2, outm, outp, outt,int,Eint,excess)
	cnames <- c("med","sdmed1","sdmed0","mean","sdmean","below","above","outliers",
		"inliers","Einliers","ExcessIns")
	dimnames(tab) <- list(paste("Day",1:24),cnames)
	}

Description

Plot an object generated by KhmaladzeTest

Usage

## S3 method for class 'KhmaladzeTest'
plot(x, ...)

Arguments

See Also

KhmaladzeTest

Description

Function to plot quantile regression process.

Usage

## S3 method for class 'rq.process'
plot(x, nrow=3, ncol=2, ...)

Arguments

Author(s)

Roger Koenker [email protected]

See Also

rq

Description

A sequence of coefficient estimates for quantile regressions with varying tau parameters is visualized.

Usage

## S3 method for class 'rqs'
plot(x, parm = NULL, ols = TRUE,
  mfrow = NULL, mar = NULL, ylim = NULL, main = NULL, col = 1:2, lty = 1:2,
  cex = 0.5, pch = 20, type = "b", xlab = "", ylab = "", ...)

Arguments

Details

The plot method for "rqs" objects visualizes the coefficients only, confidence bands can be added by using the plot method for the associated "summary.rqs" object.

Value

A matrix with all coefficients visualized is returned invisibly.

See Also

rq, plot.summary.rqs

Examples

## fit Engel models (in levels) for tau = 0.1, ..., 0.9
data("engel")
fm <- rq(foodexp ~ income, data = engel, tau = 1:9/10)

## visualizations
plot(fm)
plot(fm, parm = 2, mar = c(5.1, 4.1, 2.1, 2.1), main = "", xlab = "tau", 
  ylab = "income coefficient", cex = 1, pch = 19)

Description

Takes a fitted rqss object produced by rqss() and plots the component smooth functions that make up the ANOVA decomposition. Since the components "omit the intercept" the estimated intercept is added back in – this facilitates the comparison of quantile fits particularly. For models with a partial linear component or several qss components it may be preferable to plot the output of predict.rqss. Note that these functions are intended to plot rqss objects only, attempting to plot summary.rqss objects just generates a warning message.

Usage

## S3 method for class 'rqss'
plot(x, rug = TRUE, jit = TRUE, bands = NULL, coverage = 0.95,
	add = FALSE, shade = TRUE, select = NULL, pages = 0, titles = NULL, 
	bcol = NULL, ...)
## S3 method for class 'qss1'
plot(x, rug = TRUE, jit = TRUE, add = FALSE, ...)
## S3 method for class 'qts1'
plot(x, rug = TRUE, jit = TRUE, add = FALSE, ...)
## S3 method for class 'qss2'
plot(x, render = "contour", ncol = 100, zcol = NULL, ...)
## S3 method for class 'summary.rqss'
plot(x, ...)

Arguments

Details

For univariate qss components with Dorder = 0 the fitted function is piecewise constant, not piecewise linear. In this case the constraints are limited to increasing, decreasing or none. If bands == "uniform" then the bands are uniform bands based on the Hotelling (1939) tube approach. See also Naiman (1986), Johansen and Johnstone (1990), Sun and Loader (1994), and Krivobokova, Kneib, and Claeskens (2009), in particular the computation of the "tube length" is based on the last of these references. If bands is non null, and not "uniform" then pointwise bands are returned. Since bands for bivariate components are not (yet) supported, if requested such components will be returned as NULL.

Value

The function produces plots for the ANOVA components as a side effect. For "qss1" the "add = TRUE" can be used to overplot the fit on a scatterplot. When there are multiple pages required "par(ask = TRUE)" is turned on so that the plots may be examined sequentially. If bands != NULL then a list with three components for each qss component is returned (invisibly):

Author(s)

Roger Koenker

References

[1] Hotelling, H. (1939): “Tubes and Spheres in $n$-spaces, and a class of statistical problems,” Am J. Math, 61, 440–460.

[2] Johansen, S., and I.M. Johnstone (1990): “Hotelling's Theorem on the Volume of Tubes: Some Illustrations in Simultaneous Inference and Data Analysis,” The Annals of Statistics, 18, 652–684.

[3] Naiman, D. (1986) Conservative confidence bands in curvilinear regression, The Annals of Statistics, 14, 896–906.

[4] Sun, J. and C.R. Loader, (1994) Simultaneous confidence bands for linear regression and smoothing, The Annals of Statistics, 22, 1328–1345.

[5] Krivobokova, T., T. Kneib, and G. Claeskens (2009) Simultaneous Confidence Bands for Penalized Spline Estimators, preprint.

[6] Koenker, R. (2010) Additive Models for Quantile Regression: Model Selection and Confidence Bandaids, preprint.

See Also

rqss

Examples

n <- 200
x <- sort(rchisq(n,4))
z <- x + rnorm(n)
y <- log(x)+ .1*(log(x))^2 + log(x)*rnorm(n)/4 + z
plot(x,y-z)
fN <- rqss(y~qss(x,constraint="N")+z)
plot(fN)
fI <- rqss(y~qss(x,constraint="I")+z)
plot(fI,  col="blue")
fCI <- rqss(y~qss(x,constraint="CI")+z)
plot(fCI, col="red")

## A bivariate example
if(requireNamespace("interp")){
if(requireNamespace("interp")){
data(CobarOre)
fCO <- rqss(z~qss(cbind(x,y),lambda=.08), data = CobarOre)
plot(fCO)
}}

Description

A sequence of coefficient estimates for quantile regressions with varying tau parameters is visualized along with associated confidence bands.

Usage

## S3 method for class 'summary.rqs'
plot(x, parm = NULL, level = 0.9, ols = TRUE,
  mfrow = NULL, mar = NULL, ylim = NULL, main = NULL,
  col = gray(c(0, 0.75)), border = NULL, lcol = 2, lty = 1:2,
  cex = 0.5, pch = 20, type = "b", xlab = "", ylab = "", ...)

Arguments

Details

The plot method for "summary.rqs" objects visualizes the coefficients along with their confidence bands. The bands can be omitted by using the plot method for "rqs" objects directly.

Value

A list with components z, an array with all coefficients visualized (and associated confidence bands), and Ylim, a 2 by p matrix containing the y plotting limits. The latter component may be useful for establishing a common scale for two or more similar plots. The list is returned invisibly.

See Also

rq, plot.rqs

Examples

## fit Engel models (in levels) for tau = 0.1, ..., 0.9
data("engel")
fm <- rq(foodexp ~ income, data = engel, tau = 1:9/10)
sfm <- summary(fm)

## visualizations
plot(sfm)
plot(sfm, parm = 2, mar = c(5.1, 4.1, 2.1, 2.1), main = "", xlab = "tau", 
  ylab = "income coefficient", cex = 1, pch = 19)

Description

Prediction based on fitted quantile regression model

Usage

## S3 method for class 'rq'
predict(object, newdata, type = "none", interval = c("none", "confidence"), 
	level = .95, na.action = na.pass, ...)
## S3 method for class 'rqs'
predict(object, newdata, type = "Qhat", stepfun = FALSE, na.action = na.pass, ...)
## S3 method for class 'rq.process'
predict(object, newdata, type = "Qhat", stepfun = FALSE, na.action = na.pass, ...)

Arguments

Details

Produces predicted values, obtained by evaluating the quantile regression function in the frame 'newdata' (which defaults to 'model.frame(object)'. These predictions purport to estimate the conditional quantile function of the response variable of the fitted model evaluated at the covariate values specified in "newdata" and the quantile(s) specified by the "tau" argument. Several methods are provided to compute confidence intervals for these predictions.

Value

A vector or matrix of predictions, depending upon the setting of 'interval'. In the case that there are multiple taus in object when object is of class 'rqs' setting 'stepfun = TRUE' will produce a stepfun object or a list of stepfun objects. The function rearrange can be used to monotonize these step-functions, if desired.

Author(s)

R. Koenker

References

Zhou, Kenneth Q. and Portnoy, Stephen L. (1998) Statistical inference on heteroscedastic models based on regression quantiles Journal of Nonparametric Statistics, 9, 239-260

See Also

rq rearrange

Examples

data(airquality)
airq <- airquality[143:145,]
f <- rq(Ozone ~ ., data=airquality)
predict(f,newdata=airq)
f <- rq(Ozone ~ ., data=airquality, tau=1:19/20)
fp <- predict(f, newdata=airq, stepfun = TRUE)
fpr <- rearrange(fp)
plot(fp[[2]],main = "Conditional Ozone Quantile Prediction")
lines(fpr[[2]], col="red")
legend(.2,20,c("raw","cooked"),lty = c(1,1),col=c("black","red"))
fp <- predict(f, newdata=airq, type = "Fhat", stepfun = TRUE)
fpr <- rearrange(fp)
plot(fp[[2]],main = "Conditional Ozone Distribution Prediction")
lines(fpr[[2]], col="red")
legend(20,.4,c("raw","cooked"),lty = c(1,1),col=c("black","red"))

Description

Additive models for nonparametric quantile regression using total variation penalty methods can be fit with the rqss function. Univarariate and bivariate components can be predicted using these functions.

Usage

## S3 method for class 'rqss'
predict(object, newdata, interval = "none", level = 0.95, ...)
## S3 method for class 'qss1'
predict(object, newdata, ...)
## S3 method for class 'qss2'
predict(object, newdata, ...)

Arguments

Details

For both univariate and bivariate prediction linear interpolation is done. In the bivariate case, this involves computing barycentric coordinates of the new points relative to their enclosing triangles. It may be of interest to plot individual components of fitted rqss models: this is usually best done by fixing the values of other covariates at reference values typical of the sample data and predicting the response at varying values of one qss term at a time. Direct use of the predict.qss1 and predict.qss2 functions is discouraged since it usually corresponds to predicted values at absurd reference values of the other covariates, i.e. zero.

Value

A vector of predictions, or in the case that interval = "confidence") a matrix whose first column is the vector of predictions and whose second and third columns are the lower and upper confidence limits for each prediction.

Author(s)

R. Koenker

See Also

rqss

Examples

n <- 200
lam <- 2
x <- sort(rchisq(n,4))
z <- exp(rnorm(n)) + x
y <- log(x)+ .1*(log(x))^2 + z/4 +  log(x)*rnorm(n)/4
plot(x,y - z/4 + mean(z)/4)
Ifit <- rqss(y ~ qss(x,constraint="I") + z)
sfit <- rqss(y ~ qss(x,lambda = lam) + z)
xz <- data.frame(z = mean(z),
                 x = seq(min(x)+.01,max(x)-.01,by=.25))
lines(xz[["x"]], predict(Ifit, xz), col=2)
lines(xz[["x"]], predict(sfit, xz), col=3)
legend(10,2,c("Increasing","Smooth"),lty = 1, col = c(2,3))
title("Predicted Median Response at Mean Value of z")

## Bivariate example -- loads pkg "interp"
if(requireNamespace("interp")){
if(requireNamespace("interp")){
data(CobarOre)
fit <- rqss(z ~ qss(cbind(x,y), lambda=.08),
            data= CobarOre)
plot(fit, col="grey",
     main = "CobarOre data -- rqss(z ~ qss(cbind(x,y)))")
T <- with(CobarOre, interp::tri.mesh(x, y))
set.seed(77)
ndum <- 100
xd <- with(CobarOre, runif(ndum, min(x), max(x)))
yd <- with(CobarOre, runif(ndum, min(y), max(y)))
table(s <- interp::in.convex.hull(T, xd, yd))
pred <- predict(fit, data.frame(x = xd[s], y = yd[s]))
contour(interp::interp(xd[s],yd[s], pred),
        col="red", add = TRUE)
}}

Description

Print an object generated by KhmaladzeTest

Usage

## S3 method for class 'KhmaladzeTest'
print(x, ...)

Arguments

See Also

KhmaladzeTest

Description

Print an object generated by rq

Usage

## S3 method for class 'rq'
print(x, ...)
## S3 method for class 'rqs'
print(x, ...)

Arguments

See Also

rq

Description

Print summary of quantile regression object

Usage

## S3 method for class 'summary.rq'
print(x, digits=max(5, .Options$digits - 2), ...)
## S3 method for class 'summary.rqs'
print(x, ...)

Arguments

See Also

summary.rq

Description

The function q489 computes a single sample quantile using a fortran implementation of the Floyd and Rivest (1975) algorithm. In contrast to the more elaborate function kuantile that uses the Kiweil (2005) implementation it does not attempt to replicate the nine varieties of quantiles as documented in the base function. quantile

Usage

q489(x, tau = .5)

Arguments

Details

This is a direct translation of the Algol 68 implementation of Floyd and Rivest (1975), implemented in Ratfor. For the median, average case behavior requires $1.5 n + O((n log n)^{1/2})$ comparisons. In preliminary experiments it seems to be somewhat faster in large samples than the implementation kuantile of Kiwiel (2005). See Knuth (1998) for further details. No provision is made for non-uniqueness of the quantile. so, when $\tau n$ is an integer there may be some discrepancy.

Value

A scalar quantile of the same length as the vector p.

Author(s)

R.W.Floyd and R.L.Rivest, R implementation: Roger Koenker

References

R.W. Floyd and R.L. Rivest: "Algorithm 489: The Algorithm SELECT—for Finding the $i$th Smallest of $n$ Elements", Comm. ACM 18, 3 (1975) 173,

K.C. Kiwiel: On Floyd and Rivest's SELECT Algorithm, Theoretical Computer Sci. 347 (2005) 214-238.

D. Knuth, The Art of Computer Programming, Volume 3, Sorting and Searching, 2nd Ed., (1998), Addison-Wesley.

See Also

quantile

Examples

medx <- q489(rnorm(1001))

Description

This function solves a weighted quantile regression problem to find the optimal portfolio weights minimizing a Choquet risk criterion described in Bassett, Koenker, and Kordas (2002).

Usage

qrisk(x, alpha = c(0.1, 0.3), w = c(0.7, 0.3), mu = 0.07, 
      R = NULL, r = NULL, lambda = 10000)

Arguments

Details

The function calls rq.fit.hogg which in turn calls the constrained Frisch Newton algorithm. The constraints Rb=r are intended to apply only to the slope parameters, not the intercept parameters. The user is completely responsible to specify constraints that are consistent, ie that have at least one feasible point. See examples for imposing non-negative portfolio weights.

Value

Author(s)

R. Koenker

References

http://www.econ.uiuc.edu/~roger/research/risk/risk.html

Bassett, G., R. Koenker, G Kordas, (2004) Pessimistic Portfolio Allocation and Choquet Expected Utility, J. of Financial Econometrics, 2, 477-492.

See Also

rq.fit.hogg, srisk

Examples

#Fig 1:  ... of Choquet paper
        mu1 <- .05; sig1 <- .02; mu2 <- .09; sig2 <- .07
        x <- -10:40/100
        u <- seq(min(c(x)),max(c(x)),length=100)
        f1 <- dnorm(u,mu1,sig1)
        F1 <- pnorm(u,mu1,sig1)
        f2 <- dchisq(3-sqrt(6)*(u-mu1)/sig1,3)*(sqrt(6)/sig1)
        F2 <- pchisq(3-sqrt(6)*(u-mu1)/sig1,3)
        f3 <- dnorm(u,mu2,sig2)
        F3 <- pnorm(u,mu2,sig2)
        f4 <- dchisq(3+sqrt(6)*(u-mu2)/sig2,3)*(sqrt(6)/sig2)
        F4 <- pchisq(3+sqrt(6)*(u-mu2)/sig2,3)
        plot(rep(u,4),c(f1,f2,f3,f4),type="n",xlab="return",ylab="density")
        lines(u,f1,lty=1,col="blue")
        lines(u,f2,lty=2,col="red")
        lines(u,f3,lty=3,col="green")
        lines(u,f4,lty=4,col="brown")
        legend(.25,25,paste("Asset ",1:4),lty=1:4,col=c("blue","red","green","brown"))
#Now generate random sample of returns from these four densities.
n <- 1000
if(TRUE){ #generate a new returns sample if TRUE
	x1 <- rnorm(n)
	x1 <- (x1-mean(x1))/sqrt(var(x1))
	x1 <- x1*sig1 + mu1
	x2 <- -rchisq(n,3)
	x2 <- (x2-mean(x2))/sqrt(var(x2))
	x2 <- x2*sig1 +mu1
	x3 <- rnorm(n)
	x3 <- (x3-mean(x3))/sqrt(var(x3))
	x3 <- x3*sig2 +mu2
	x4 <- rchisq(n,3)
	x4 <- (x4-mean(x4))/sqrt(var(x4))
	x4 <- x4*sig2 +mu2
	}
library(quantreg)
x <- cbind(x1,x2,x3,x4)
qfit <- qrisk(x)
sfit <- srisk(x)
# Try new distortion function
qfit1 <- qrisk(x,alpha = c(.05,.1), w = c(.9,.1),mu = 0.09)
# Constrain portfolio weights to be non-negative
qfit2 <- qrisk(x,alpha = c(.05,.1), w = c(.9,.1),mu = 0.09,
	       R = rbind(rep(-1,3), diag(3)), r = c(-1, rep(0,3)))

Description

In the formula specification of rqss nonparametric terms are specified with qss. Both univariate and bivariate specifications are possible, and qualitative constraints may also be specified for the qss terms.

Usage

qss(x, constraint = "N", lambda = 1, ndum = 0, dummies = NULL, 
    Dorder = 1, w = rep(1, length(x)))

Arguments

Details

The various pieces returned are stored in sparse matrix.csr form. See rqss for details on how they are assembled. To preserve the sparsity of the design matrix the first column of each qss term is dropped. This differs from the usual convention that would have forced qss terms to have mean zero. This convention has implications for prediction that need to be recognized. The penalty components for qss terms are based on total variation penalization of the first derivative (and gradient, for bivariate x) as described in the references appearing in the help for rqss. When Dorder = 0, fitting is like the taut string methods of Davies (2014), except for the fact that fidelity is quantilesque rather than quadratic, and that no provision is made for automatic selection of the smoothing parameter.

For the bivariate case, package interp (and for plotting also interp) are required (automatically, by the R code).

Value

Author(s)

Roger Koenker

References

Davies, Laurie (2014) Data Analysis and Approximate Models, CRC Press.

See Also

rqss

Description

Computes quantile treatment effects comparable to those of crq model from a coxph object.

Usage

QTECox(x, smooth = TRUE)

Arguments

Details

Estimates of the Cox QTE, $\frac{dQ(t|x)}{dx_{j}}$ at $x=\bar{x}$, can be expressed as a function of t as follows:

$\frac{dQ(t|x)}{dx_{j}}=\frac{dt}{dx_{j}}\frac{dQ(t|x)}{dt}$

The Cox survival function, $S(y|x)=\exp \{-H_{0}(y)\exp (b^{\prime }x)\}$

$\frac{dS(y|x)}{dx_{j}}=S(y|x)log \{S(y|x)\}b_{j}$

where $\frac{dQ(t|x)}{dx_{j}}$ can be estimated by $\frac{\Delta (t)}{\Delta (S)} (1-t)$ where $S$ and $t$ denote the surv and time components of the survfit object. Note that since $t=1-S(y|x)$, the above is the value corresponding to the argument $(1-t)$; and furthermore

$\frac{dt}{dx_{j}}=-\frac{dS(y|x)}{dx_{j}}=-(1-t) log (1-t)b_{j}$

Thus the QTE at the mean of x's is:

$(1-S)= \frac{\Delta (t)}{\Delta (S)}S ~log (S)b_{j}$

Since $\Delta S$ is negative and $log (S)$ is also negative this has the same sign as $b_{j}$ The crq model fits the usual AFT form Surv(log(Time),Status), then

$\frac{d log (Q(t|x))}{dx_{j}}=\frac{dQ(t|x)}{dx_{j}}/ Q(t|x)$

This is the matrix form returned.

Value

Author(s)

Roger Koenker Stephen Portnoy & Tereza Neocleous

References

Koenker, R. and Geling, O. (2001). Reappraising Medfly longevity: a quantile regression survival analysis, J. Amer. Statist. Assoc., 96, 458-468

See Also

crq

Description

Function to compute ranks from the dual (regression rankscore) process.

Usage

ranks(v, score="wilcoxon", tau=0.5, trim = NULL)

Arguments

Details

See GJKP(1993) for further details.

Value

The function returns two components. One is the ranks, the other is a scale factor which is the $L_2$ norm of the score function. All score functions should be normalized to have mean zero.