R : Copyright 2005, The R Foundation for Statistical Computing Version 2.1.1 (2005-06-20), ISBN 3-900051-07-0 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for a HTML browser interface to help. Type 'q()' to quit R. > ### *
> ### > attach(NULL, name = "CheckExEnv") > assign(".CheckExEnv", as.environment(2), pos = length(search())) # base > ## add some hooks to label plot pages for base and grid graphics > setHook("plot.new", ".newplot.hook") > setHook("persp", ".newplot.hook") > setHook("grid.newpage", ".gridplot.hook") > > assign("cleanEx", + function(env = .GlobalEnv) { + rm(list = ls(envir = env, all.names = TRUE), envir = env) + RNGkind("default", "default") + set.seed(1) + options(warn = 1) + delayedAssign("T", stop("T used instead of TRUE"), + assign.env = .CheckExEnv) + delayedAssign("F", stop("F used instead of FALSE"), + assign.env = .CheckExEnv) + sch <- search() + newitems <- sch[! sch %in% .oldSearch] + for(item in rev(newitems)) + eval(substitute(detach(item), list(item=item))) + missitems <- .oldSearch[! .oldSearch %in% sch] + if(length(missitems)) + warning("items ", paste(missitems, collapse=", "), + " have been removed from the search path") + }, + env = .CheckExEnv) > assign("..nameEx", "__{must remake R-ex/*.R}__", env = .CheckExEnv) # for now > assign("ptime", proc.time(), env = .CheckExEnv) > grDevices::postscript("quantreg-Examples.ps") > assign("par.postscript", graphics::par(no.readonly = TRUE), env = .CheckExEnv) > options(contrasts = c(unordered = "contr.treatment", ordered = "contr.poly")) > options(warn = 1) > library('quantreg') quantreg package loaded quantreg now includes the nlrq and nprq packages > > assign(".oldSearch", search(), env = .CheckExEnv) > assign(".oldNS", loadedNamespaces(), env = .CheckExEnv) > cleanEx(); ..nameEx <- "CobarOre" > > ### * CobarOre > > flush(stderr()); flush(stdout()) > > ### Name: CobarOre > ### Title: Cobar Ore data > ### Aliases: CobarOre > ### Keywords: datasets > > ### ** Examples > > data(CobarOre) > plot(CobarOre) > > > > cleanEx(); ..nameEx <- "anova.rq" > > ### * anova.rq > > flush(stderr()); flush(stdout()) > > ### Name: anova.rq > ### Title: Anova function for quantile regression fits > ### Aliases: anova.rq anova.rqlist print.anova.rq rq.test.rank > ### Keywords: htest regression robust > > ### ** 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) Quantile Regression Analysis of Variance Table Model 1: y.net ~ lgdp2 + fse2 + gedy2 + Iy2 + gcony2 Model 2: y.net ~ lgdp2 + fse2 + gedy2 Df Resid Df F value Pr(>F) 1 2 155 18.879 4.596e-08 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > anova(fit1,fit2,fit3) Warning in summary.rq(x, se = "nid") : 1 non-positive fis Quantile Regression Analysis of Variance Table Model: y.net ~ lgdp2 + fse2 + gedy2 + Iy2 + gcony2 Joint Test of Equality of Slopes: tau in { 0.5 0.75 0.25 } Df Resid Df F value Pr(>F) 1 10 473 1.8039 0.05751 . --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > anova(fit1,fit2,fit3,joint=FALSE) Warning in summary.rq(x, se = "nid") : 1 non-positive fis Quantile Regression Analysis of Variance Table Model: y.net ~ lgdp2 + fse2 + gedy2 + Iy2 + gcony2 Tests of Equality of Distinct Slopes: tau in { 0.5 0.75 0.25 } Df Resid Df F value Pr(>F) lgdp2 2 481 1.0656 0.3453 fse2 2 481 2.6399 0.0724 . gedy2 2 481 0.7862 0.4561 Iy2 2 481 0.0447 0.9563 gcony2 2 481 0.0653 0.9368 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > > > > cleanEx(); ..nameEx <- "boot.rq" > > ### * boot.rq > > flush(stderr()); flush(stdout()) > > ### Name: boot.rq > ### Title: Bootstrapping Quantile Regression > ### Aliases: boot.rq boot.rq.xy boot.rq.pwy boot.rq.mcmb > ### Keywords: regression > > ### ** Examples > > y <- rnorm(50) > x <- matrix(rnorm(100),50) > fit <- rq(y~x,tau = .4) > summary(fit,se = "boot", bsmethod= "xy") Call: rq(formula = y ~ x, tau = 0.4) tau: [1] 0.4 Coefficients: Value Std. Error t value Pr(>|t|) (Intercept) -0.07369 0.16975 -0.43414 0.66617 x1 -0.18160 0.16929 -1.07270 0.28888 x2 0.09891 0.16081 0.61508 0.54147 > summary(fit,se = "boot", bsmethod= "pwy") Call: rq(formula = y ~ x, tau = 0.4) tau: [1] 0.4 Coefficients: Value Std. Error t value Pr(>|t|) (Intercept) -0.07369 0.17008 -0.43330 0.66678 x1 -0.18160 0.17942 -1.01217 0.31664 x2 0.09891 0.20797 0.47560 0.63656 > #summary(fit,se = "boot", bsmethod= "mcmb") > > > > cleanEx(); ..nameEx <- "khmaladze.test" > > ### * khmaladze.test > > flush(stderr()); flush(stdout()) > > ### Name: khmaladze.test > ### Title: Tests of Location and Location Scale Hypothesis for Linear > ### Models > ### Aliases: khmaladze.test > ### Keywords: htest > > ### ** Examples > > data(barro) > fit <- rqProcess( y.net ~ lgdp2 + fse2 + gedy2 + Iy2 + gcony2, + data = barro, taus = seq(.1,.9,by = .05)) taus: 0.1 Warning in summary.rq(rq(formula, tau = taus[i], method = "fn", data = data), : 1 non-positive fis 0.15 Warning in summary.rq(rq(formula, tau = taus[i], method = "fn", data = data), : 3 non-positive fis 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 Warning in summary.rq(rq(formula, tau = taus[i], method = "fn", data = data), : 1 non-positive fis 0.65 Warning in summary.rq(rq(formula, tau = taus[i], method = "fn", data = data), : 1 non-positive fis 0.7 0.75 0.8 0.85 0.9 > khmaladze.test(fit, nullH = "location") $nullH [1] "location" $Tn [1] 5.221146 $THn (Intercept) lgdp2 fse2 gedy2 Iy2 gcony2 0.05043884 0.69436390 1.21528613 0.78140013 0.70499000 1.36268796 attr(,"class") [1] "khmaladze" > > > > cleanEx(); ..nameEx <- "latex.summary.rqs" > > ### * latex.summary.rqs > > flush(stderr()); flush(stdout()) > > ### Name: latex.summary.rqs > ### Title: Make a latex table from a table of rq results > ### Aliases: latex.summary.rqs > ### Keywords: IO > > ### ** Examples > > > > > cleanEx(); ..nameEx <- "latex.table" > > ### * latex.table > > flush(stderr()); flush(stdout()) > > ### Name: latex.table > ### Title: Writes a latex formatted table to a file > ### Aliases: latex.table > ### Keywords: utilities > > ### ** Examples > > > > > cleanEx(); ..nameEx <- "lprq" > > ### * lprq > > flush(stderr()); flush(stdout()) > > ### Name: lprq > ### Title: locally polynomial quantile regression > ### Aliases: lprq > ### Keywords: smooth robust > > ### ** Examples > > require(MASS) Loading required package: MASS [1] TRUE > 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)) > > > > cleanEx(); ..nameEx <- "nlrq" > > ### * nlrq > > flush(stderr()); flush(stdout()) > > ### Name: nlrq > ### Title: Function to compute nonlinear quantile regression estimates > ### Aliases: nlrq nlrqModel print.nlrq summary.nlrq deviance.nlrq > ### formula.nlrq coef.nlrq fitted.nlrq predict.nlrq print.summary.nlrq > ### tau.nlrq > ### Keywords: models regression nonlinear robust > > ### ** 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 Nonlinear regression model model: y ~ SSlogis(x, Asym, mid, scal) data: Dat Asym mid scal 9.968029 11.947209 1.962115 residual sum-of-squares: 241.7909 > 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) 109.0590 : 9.968029 11.947209 1.962115 final value 108.942726 converged lambda = 1 108.9427 : 9.958648 11.943273 1.967144 iter 0 value 108.490958 final value 108.490939 converged lambda = 0.9750981 108.4909 : 9.949430 11.987472 1.998607 iter 0 value 108.471416 final value 108.471416 converged lambda = 0.99993 108.4714 : 9.941630 11.990771 1.993440 final value 108.471243 converged lambda = 1 108.4712 : 9.941008 11.990550 1.992921 iter 0 value 108.471202 final value 108.470935 converged lambda = 0.8621249 108.4709 : 9.942734 11.992773 1.993209 iter 0 value 108.470923 final value 108.470923 converged lambda = 0.9999613 108.4709 : 9.942629 11.992728 1.993136 final value 108.470919 converged lambda = 1 108.4709 : 9.942644 11.992737 1.993144 final value 108.470919 converged lambda = 1 108.4709 : 9.942644 11.992737 1.993144 final value 108.470919 converged lambda = 1 108.4709 : 9.942644 11.992737 1.993144 > 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) 108.6656 : 9.968029 11.947209 1.962115 final value 89.108243 converged lambda = 1 89.10824 : 9.432250 11.803924 1.923472 final value 85.688895 converged lambda = 1 85.6889 : 9.183598 11.794244 1.929699 iter 0 value 85.477812 final value 85.473712 converged lambda = 0.6405076 85.47371 : 9.212527 11.844090 1.938003 final value 85.447786 converged lambda = 1 85.44779 : 9.234097 11.863975 1.949241 final value 85.446407 converged lambda = 1 85.4464 : 9.242009 11.866644 1.954192 final value 85.445691 converged lambda = 1 85.44569 : 9.234247 11.864554 1.952338 final value 85.444920 converged lambda = 1 85.44492 : 9.232975 11.863979 1.953587 iter 0 value 85.448100 final value 85.443854 converged lambda = 0.3632370 85.44385 : 9.233661 11.864280 1.957197 iter 0 value 85.445154 final value 85.443668 converged lambda = 0.8495473 85.44367 : 9.233453 11.860020 1.957831 iter 0 value 85.444111 final value 85.443667 converged lambda = 0.008522615 85.44367 : 9.233449 11.860007 1.957814 final value 85.443584 converged lambda = 1 85.44358 : 9.232996 11.859020 1.955928 iter 0 value 85.443586 final value 85.443586 converged lambda = 0.9999957 85.44359 : 9.232995 11.859024 1.955916 > 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) 109.4525 : 9.968029 11.947209 1.962115 final value 89.561438 converged lambda = 1 89.56144 : 10.64021 12.13202 2.02044 final value 87.302043 converged lambda = 1 87.30204 : 10.652294 11.966018 1.958371 final value 87.200715 converged lambda = 1 87.20072 : 10.666754 11.953497 1.962447 iter 0 value 87.134039 final value 87.131462 stopped after 4 iterations lambda = 0.865945 87.13146 : 10.639094 11.949236 1.971242 iter 0 value 87.129208 final value 87.125796 stopped after 5 iterations lambda = 0.6273177 87.1258 : 10.647783 11.962634 1.975851 iter 0 value 87.123519 final value 87.122717 converged lambda = 0.804119 87.12272 : 10.647957 11.963190 1.973657 final value 87.121592 converged lambda = 1 87.1216 : 10.649877 11.962363 1.973516 final value 87.121427 converged lambda = 1 87.12143 : 10.649051 11.961685 1.973086 iter 0 value 87.121787 final value 87.121355 converged lambda = 0.5468911 87.12135 : 10.648209 11.961208 1.972643 iter 0 value 87.121400 final value 87.121400 converged lambda = 0.9999045 87.1214 : 10.648073 11.961122 1.972568 > 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) 108.3115 : 9.968029 11.947209 1.962115 final value 62.166612 converged lambda = 1 62.16661 : 9.432250 11.803924 1.923472 final value 16.887325 converged lambda = 1 16.88732 : 8.006640 11.718631 1.979243 iter 0 value 15.823507 final value 15.823276 converged lambda = 0.7133884 15.82328 : 8.135460 12.048708 1.987995 iter 0 value 15.732947 final value 15.732737 converged lambda = 0.7726633 15.73274 : 8.042059 12.019441 1.994386 iter 0 value 15.748295 final value 15.732737 converged lambda = 0 15.73274 : 8.042059 12.019441 1.994386 > 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) 109.8066 : 9.968029 11.947209 1.962115 final value 56.575819 converged lambda = 1 56.57582 : 10.672415 12.148657 2.027285 final value 20.551829 converged lambda = 1 20.55183 : 11.923558 12.366710 2.121476 final value 17.268734 converged lambda = 1 17.26873 : 12.266850 12.051876 2.060768 iter 0 value 17.229809 final value 17.194623 converged lambda = 0.5512476 17.19462 : 12.176373 12.020546 2.003537 iter 0 value 17.177793 final value 17.175845 stopped after 3 iterations lambda = 0.9001389 17.17585 : 12.180837 12.005129 2.019783 iter 0 value 17.187685 final value 17.175761 converged lambda = 0.1504769 17.17576 : 12.177202 12.003960 2.011709 final value 17.175612 converged lambda = 1 17.17561 : 12.181536 12.005336 2.018940 final value 17.175603 converged lambda = 1 17.17560 : 12.181679 12.005403 2.019175 final value 17.175518 converged lambda = 1 17.17552 : 12.179538 12.004686 2.014531 iter 0 value 17.243132 final value 17.175518 converged lambda = 0 17.17552 : 12.179538 12.004686 2.014531 > 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) > > > > cleanEx(); ..nameEx <- "plot" > > ### * plot > > flush(stderr()); flush(stdout()) > > ### Name: plot > ### Title: Default Ploting Method for rqss() > ### Aliases: plot.rqss plot.qss1 plot.qss2 > ### Keywords: regression smooth iplot > > ### ** 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) > fit <- rqss(y~qss(x,constraint="N")+z) Loading required package: SparseM [1] "SparseM library loaded" > lines(x[-1],fit$coef[1]+fit$coef[-(1:2)]) > fit <- rqss(y~qss(x,constraint="I")+z) > lines(x[-1],fit$coef[1]+fit$coef[-(1:2)],col="blue") > fit <- rqss(y~qss(x,constraint="CI")+z) > lines(x[-1],fit$coef[1]+fit$coef[-(1:2)],col="red") > #Cleanup > rm(list=ls()) > #A bivariate example > data(CobarOre) > attach(CobarOre) > fit <- rqss(z~qss(cbind(x,y),lambda=.08)) Loading required package: tripack > plot(fit) Loading required package: akima > > > > cleanEx(); ..nameEx <- "predict.qss" > > ### * predict.qss > > flush(stderr()); flush(stdout()) > > ### Name: predict.qss > ### Title: Predict based on nonparametric quantile regression smoothing > ### spline component > ### Aliases: predict.qss predict.qss1 predict.qss2 > ### Keywords: regression smooth robust > > ### ** Examples > > > > cleanEx(); ..nameEx <- "predict.rq" > > ### * predict.rq > > flush(stderr()); flush(stdout()) > > ### Name: predict.rq > ### Title: Quantile Regression Prediction > ### Aliases: predict.rq predict.rqs > ### Keywords: regression > > ### ** Examples > > data(airquality) > airq <- airquality[143:150,] > f <- rq(Ozone ~ ., data=airquality) > predict(f,newdata=airq) 143 144 145 146 147 148 149 48.963525 2.721109 17.638599 39.170443 12.951419 -17.991994 31.854806 150 25.103796 > f <- rq(Ozone ~ ., data=airquality,tau=1:3/4) > predict(f,newdata=airq) tau= 0.25 tau= 0.50 tau= 0.75 143 43.907687 48.963525 59.06222 144 1.169585 2.721109 10.02194 145 11.116311 17.638599 23.11926 146 32.757548 39.170443 49.14094 147 7.000000 12.951419 17.77711 148 -23.618139 -17.991994 -10.43453 149 27.779496 31.854806 33.57666 150 18.785426 25.103796 34.13604 > > > > cleanEx(); ..nameEx <- "ranks" > > ### * ranks > > flush(stderr()); flush(stdout()) > > ### Name: ranks > ### Title: Quantile Regression Ranks > ### Aliases: ranks > ### Keywords: regression > > ### ** Examples > > data(stackloss) > ranks(rq(stack.loss ~ stack.x, tau=-1)) $ranks [1] 0.18302994 -0.01602686 0.35483163 0.45653019 -0.19694479 -0.34096965 [7] -0.23121754 -0.01788246 -0.43795303 0.11192103 0.23876556 0.18554420 [13] -0.29336169 -0.17945242 0.36809323 0.02610216 -0.22538240 -0.01259556 [19] 0.16445517 0.29778908 -0.43527579 $A2 [1] 0.08333333 > > > > cleanEx(); ..nameEx <- "rq" > > ### * rq > > flush(stderr()); flush(stdout()) > > ### Name: rq > ### Title: Quantile Regression > ### Aliases: rq > ### Keywords: regression > > ### ** Examples > > data(stackloss) > rq(stack.loss ~ stack.x,.5) #median (l1) regression fit for the stackloss data. Call: rq(formula = stack.loss ~ stack.x, tau = 0.5) Coefficients: (Intercept) stack.xAir.Flow stack.xWater.Temp stack.xAcid.Conc. -39.68985507 0.83188406 0.57391304 -0.06086957 Degrees of freedom: 21 total; 17 residual > rq(stack.loss ~ stack.x,.25) #the 1st quartile, Call: rq(formula = stack.loss ~ stack.x, tau = 0.25) Coefficients: (Intercept) stack.xAir.Flow stack.xWater.Temp stack.xAcid.Conc. -3.60000e+01 5.00000e-01 1.00000e+00 -4.57967e-16 Degrees of freedom: 21 total; 17 residual > #note that 8 of the 21 points lie exactly on this plane in 4-space! > rq(stack.loss ~ stack.x, tau=-1) #this returns the full rq process $sol [,1] [,2] [,3] [,4] tau 0.00000000 0.1240939 1.300537e-01 1.303646e-01 Qbar 0.00000000 13.9936756 1.530952e+01 1.530952e+01 Obj.Fun 0.00000000 10.6056828 1.104750e+01 1.106196e+01 (Intercept) -29.01401869 -36.0781250 -3.600000e+01 -3.600000e+01 stack.xAir.Flow 0.31542056 0.3515625 5.000000e-01 5.000000e-01 stack.xWater.Temp 1.22429907 1.7500000 1.000000e+00 1.000000e+00 stack.xAcid.Conc. -0.02803738 -0.0937500 -1.443290e-15 -4.758099e-16 [,5] [,6] [,7] [,8] tau 1.494253e-01 1.607143e-01 1.945662e-01 2.231167e-01 Qbar 1.530952e+01 1.530952e+01 1.530952e+01 1.530952e+01 Obj.Fun 1.194828e+01 1.247321e+01 1.404733e+01 1.537493e+01 (Intercept) -3.600000e+01 -3.600000e+01 -3.600000e+01 -3.600000e+01 stack.xAir.Flow 5.000000e-01 5.000000e-01 5.000000e-01 5.000000e-01 stack.xWater.Temp 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 stack.xAcid.Conc. -4.758099e-16 -1.400996e-16 -8.679925e-17 -4.758099e-16 [,9] [,10] [,11] [,12] tau 2.539683e-01 0.27510618 0.3310042 0.3749882 Qbar 1.530952e+01 16.16141457 16.4441422 16.8013416 Obj.Fun 1.680952e+01 17.79243732 19.3916975 20.3889465 (Intercept) -3.600000e+01 -37.89705882 -38.5266594 -32.6377709 stack.xAir.Flow 5.000000e-01 0.75735294 0.8427639 0.8250774 stack.xWater.Temp 1.000000e+00 0.79411765 0.7257889 0.7399381 stack.xAcid.Conc. -4.758099e-16 -0.09803922 -0.1305767 -0.1857585 [,13] [,14] [,15] [,16] tau 0.3918757 0.40948814 0.48984468 0.56478766 Qbar 16.9593496 17.42450639 17.43436853 17.44517734 Obj.Fun 20.6451613 20.85393258 21.02150538 21.16226784 (Intercept) -33.2461197 -39.65188470 -39.68985507 -39.73147023 stack.xAir.Flow 0.8492239 0.83037694 0.83188406 0.83353584 stack.xWater.Temp 0.6274945 0.58093126 0.57391304 0.56622114 stack.xAcid.Conc. -0.1662971 -0.06208426 -0.06086957 -0.05953827 [,17] [,18] [,19] [,20] tau 0.5923717 0.604223291 0.619988864 0.6511309 Qbar 17.4565950 19.136255092 19.137513359 19.1484286 Obj.Fun 21.2078167 21.224545264 20.690701559 19.6353667 (Intercept) -39.6814159 -54.032448378 -54.064837905 -54.0700000 stack.xAir.Flow 0.8362832 0.873156342 0.872817955 0.8710000 stack.xWater.Temp 0.5663717 0.979351032 0.980049875 0.9840000 stack.xAcid.Conc. -0.0619469 -0.002949853 -0.002493766 -0.0020000 [,21] [,22] [,23] [,24] tau 6.897262e-01 7.621009e-01 0.76843239 0.77392070 Qbar 1.915640e+01 1.919264e+01 19.71524260 19.98904006 Obj.Fun 1.831861e+01 1.583728e+01 15.61539215 15.36281938 (Intercept) -5.418966e+01 -5.409091e+01 -54.33806147 -56.68253968 stack.xAir.Flow 8.706897e-01 8.636364e-01 0.77659574 0.78571429 stack.xWater.Temp 9.827586e-01 1.000000e+00 1.18912530 1.25396825 stack.xAcid.Conc. -5.698567e-16 -5.843653e-16 0.02364066 0.03174603 [,25] [,26] [,27] [,28] [,29] tau 0.77767777 0.8142857 0.83392070 0.9130604 1.0000000 Qbar 20.12133072 20.1607143 20.20634921 21.7007310 0.0000000 Obj.Fun 15.16831683 13.1714286 12.08414097 7.6259394 0.0000000 (Intercept) -58.54794521 -59.3750000 -58.54331865 -58.4619970 -58.4619970 stack.xAir.Flow 0.80821918 0.8062500 0.79295154 0.5245902 0.5245902 stack.xWater.Temp 1.27397260 1.2562500 1.30543319 1.8584203 1.8584203 stack.xAcid.Conc. 0.03424658 0.0500000 0.03817915 0.1073025 0.1073025 $dsol [,1] [,2] [,3] [,4] [,5] [,6] [,7] [1,] 1 1.00000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 [2,] 1 1.00000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 [3,] 1 1.00000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 [4,] 1 1.00000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 [5,] 1 1.00000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 [6,] 1 1.00000000 1.0000000 0.9541788 0.6063218 0.4598214 0.0000000 [7,] 1 1.00000000 0.9598530 1.0000000 1.0000000 1.0000000 1.0000000 [8,] 1 1.00000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 [9,] 1 0.00000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 [10,] 1 1.00000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 [11,] 1 1.00000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 [12,] 1 1.00000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 [13,] 1 1.00000000 1.0000000 1.0000000 1.0000000 0.9107143 0.6737511 [14,] 1 1.00000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 [15,] 1 1.00000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 [16,] 1 1.00000000 1.0000000 1.0000000 1.0000000 1.0000000 0.8871604 [17,] 1 0.85416063 0.8306474 0.8609229 0.9589491 1.0000000 1.0000000 [18,] 1 1.00000000 1.0000000 1.0000000 0.2967980 0.2544643 0.3531989 [19,] 1 0.49869527 0.4783715 0.4472410 1.0000000 1.0000000 1.0000000 [20,] 1 1.00000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 [21,] 1 0.04117135 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 [,8] [,9] [,10] [,11] [,12] [,13] [,14] [1,] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 [2,] 1.0000000 1.0000000 1.0000000 1.0000000 0.8253325 0.7371565 0.7253433 [3,] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 [4,] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 [5,] 1.0000000 1.0000000 1.0000000 0.0000000 0.0000000 0.0000000 0.0000000 [6,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 [7,] 0.7038528 0.3650794 0.1790656 0.2871626 0.1795718 0.1720430 0.0000000 [8,] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 [9,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 [10,] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 0.8202247 [11,] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 [12,] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 [13,] 0.3364002 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 [14,] 1.0000000 1.0000000 0.7003699 0.5080148 0.1203433 0.0000000 0.0000000 [15,] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 [16,] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 [17,] 1.0000000 0.8526077 0.3433347 0.2537334 0.0000000 0.0000000 0.0000000 [18,] 0.2742956 0.4489796 1.0000000 1.0000000 1.0000000 0.8614098 0.8551810 [19,] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 [20,] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 [21,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 [,15] [,16] [,17] [,18] [,19] [,20] [1,] 1.0000000 1.0000000 1.000000e+00 1.0000000 8.887806e-01 0.6723783 [2,] 0.6714456 0.1067666 0.000000e+00 0.0000000 0.000000e+00 0.0000000 [3,] 1.0000000 1.0000000 1.000000e+00 1.0000000 1.000000e+00 1.0000000 [4,] 1.0000000 1.0000000 1.000000e+00 1.0000000 1.000000e+00 1.0000000 [5,] 0.0000000 0.0000000 0.000000e+00 0.0000000 0.000000e+00 0.0000000 [6,] 0.0000000 0.0000000 0.000000e+00 0.0000000 0.000000e+00 0.0000000 [7,] 0.0000000 0.0000000 0.000000e+00 0.0000000 0.000000e+00 0.0000000 [8,] 0.2150538 0.2590420 1.254120e-01 0.0000000 0.000000e+00 0.0000000 [9,] 0.0000000 0.0000000 0.000000e+00 0.0000000 0.000000e+00 0.0000000 [10,] 0.0000000 0.0000000 0.000000e+00 0.0000000 3.439614e-10 0.2789582 [11,] 1.0000000 1.0000000 1.000000e+00 1.0000000 1.000000e+00 1.0000000 [12,] 1.0000000 1.0000000 8.457071e-01 0.7093874 6.948775e-01 0.3749143 [13,] 0.0000000 0.0000000 0.000000e+00 0.0000000 0.000000e+00 0.0000000 [14,] 0.0000000 0.0000000 0.000000e+00 0.0000000 0.000000e+00 0.0000000 [15,] 1.0000000 1.0000000 1.000000e+00 1.0000000 1.000000e+00 1.0000000 [16,] 1.0000000 0.0000000 5.836139e-10 0.1670500 0.000000e+00 0.0000000 [17,] 0.0000000 0.0000000 0.000000e+00 0.0000000 0.000000e+00 0.0000000 [18,] 0.8267622 0.7736505 5.890755e-01 0.4348735 3.965757e-01 0.0000000 [19,] 1.0000000 1.0000000 1.000000e+00 1.0000000 1.000000e+00 1.0000000 [20,] 1.0000000 1.0000000 1.000000e+00 1.0000000 1.000000e+00 1.0000000 [21,] 0.0000000 0.0000000 0.000000e+00 0.0000000 0.000000e+00 0.0000000 [,21] [,22] [,23] [,24] [,25] [,26] [1,] 0.4395718 5.149331e-02 0.03681582 0.0000000 0.0000000 0.0000000 [2,] 0.0000000 0.000000e+00 0.00000000 0.0000000 0.0000000 0.0000000 [3,] 1.0000000 1.000000e+00 1.00000000 1.0000000 1.0000000 0.8142857 [4,] 1.0000000 1.000000e+00 1.00000000 1.0000000 1.0000000 1.0000000 [5,] 0.0000000 0.000000e+00 0.00000000 0.0000000 0.0000000 0.0000000 [6,] 0.0000000 0.000000e+00 0.00000000 0.0000000 0.0000000 0.0000000 [7,] 0.0000000 0.000000e+00 0.00000000 0.0000000 0.0000000 0.0000000 [8,] 0.0000000 0.000000e+00 0.00000000 0.0000000 0.0000000 0.0000000 [9,] 0.0000000 0.000000e+00 0.00000000 0.0000000 0.0000000 0.0000000 [10,] 0.4703521 6.178052e-01 0.54665622 0.5090749 0.6273627 0.9000000 [11,] 1.0000000 3.265820e-01 0.00000000 0.0000000 0.0000000 0.0000000 [12,] 0.0000000 1.715086e-09 0.27944776 0.3048458 0.2106211 0.0000000 [13,] 0.0000000 0.000000e+00 0.00000000 0.0000000 0.0000000 0.0000000 [14,] 0.0000000 0.000000e+00 0.00000000 0.0000000 0.0000000 0.0000000 [15,] 1.0000000 1.000000e+00 1.00000000 0.9337445 1.0000000 1.0000000 [16,] 0.0000000 0.000000e+00 0.00000000 0.0000000 0.0000000 0.0000000 [17,] 0.0000000 0.000000e+00 0.00000000 0.0000000 0.0000000 0.0000000 [18,] 0.0000000 0.000000e+00 0.00000000 0.0000000 0.0000000 0.0000000 [19,] 0.6058266 0.000000e+00 0.00000000 0.0000000 0.0000000 0.0000000 [20,] 1.0000000 1.000000e+00 1.00000000 1.0000000 0.8307831 0.1857143 [21,] 0.0000000 0.000000e+00 0.00000000 0.0000000 0.0000000 0.0000000 [,27] [,28] [,29] [1,] 0.000000e+00 0.00000000 0 [2,] 0.000000e+00 0.00000000 0 [3,] 6.850220e-01 0.02552475 0 [4,] 1.000000e+00 1.00000000 0 [5,] 0.000000e+00 0.00000000 0 [6,] 0.000000e+00 0.00000000 0 [7,] 0.000000e+00 0.00000000 0 [8,] 0.000000e+00 0.00000000 0 [9,] 0.000000e+00 0.00000000 0 [10,] 9.057269e-01 0.27260948 0 [11,] 0.000000e+00 0.00000000 0 [12,] 2.573799e-10 0.52759782 0 [13,] 0.000000e+00 0.00000000 0 [14,] 0.000000e+00 0.00000000 0 [15,] 8.969163e-01 0.00000000 0 [16,] 0.000000e+00 0.00000000 0 [17,] 0.000000e+00 0.00000000 0 [18,] 0.000000e+00 0.00000000 0 [19,] 0.000000e+00 0.00000000 0 [20,] 0.000000e+00 0.00000000 0 [21,] 0.000000e+00 0.00000000 0 $fitted.values numeric(0) $formula stack.loss ~ stack.x $terms stack.loss ~ stack.x attr(,"variables") list(stack.loss, stack.x) attr(,"factors") stack.x stack.loss 0 stack.x 1 attr(,"term.labels") [1] "stack.x" attr(,"order") [1] 1 attr(,"intercept") [1] 1 attr(,"response") [1] 1 attr(,".Environment") attr(,"predvars") list(stack.loss, stack.x) attr(,"dataClasses") stack.loss stack.x "numeric" "nmatrix.3" $call rq(formula = stack.loss ~ stack.x, tau = -1) $tau [1] -1 $model stack.loss stack.x.Air.Flow stack.x.Water.Temp stack.x.Acid.Conc. 1 42 80 27 89 2 37 80 27 88 3 37 75 25 90 4 28 62 24 87 5 18 62 22 87 6 18 62 23 87 7 19 62 24 93 8 20 62 24 93 9 15 58 23 87 10 14 58 18 80 11 14 58 18 89 12 13 58 17 88 13 11 58 18 82 14 12 58 19 93 15 8 50 18 89 16 7 50 18 86 17 8 50 19 72 18 8 50 19 79 19 9 50 20 80 20 15 56 20 82 21 15 70 20 91 attr(,"class") [1] "rq.process" > rq(rnorm(50) ~ 1, ci=FALSE) #ordinary sample median --no rank inversion ci Warning in rq.fit.br(x, y, tau = tau, ...) : Solution may be nonunique Call: rq(formula = rnorm(50) ~ 1, ci = FALSE) Coefficients: (Intercept) 0.07456498 Degrees of freedom: 50 total; 49 residual > rq(rnorm(50) ~ 1, weights=runif(50),ci=FALSE) #weighted sample median Call: rq(formula = rnorm(50) ~ 1, weights = runif(50), ci = FALSE) Coefficients: (Intercept) 0.1532533 Degrees of freedom: 50 total; 49 residual > #plot of engel data and some rq lines see KB(1982) for references to data > data(engel) > attach(engel) > plot(x,y,xlab="Household Income",ylab="Food Expenditure",type = "n", cex=.5) > points(x,y,cex=.5,col="blue") > taus <- c(.05,.1,.25,.75,.9,.95) > xx <- seq(min(x),max(x),100) > f <- coef(rq((y)~(x),tau=taus)) > yy <- cbind(1,xx)%*%f > for(i in 1:length(taus)){ + lines(xx,yy[,i],col = "gray") + } > abline(lm(y~x),col="red",lty = 2) > abline(rq(y~x),col="blue") > legend(3000,500,c("mean (LSE) fit", "median (LAE) fit"),col = c("red","blue"),lty = c(2,1)) > #Example of plotting of coefficients and their confidence bands > plot(summary(rq(y~x,tau = 1:49/50,data=engel))) Warning in rq.fit.br(x, y, tau = tau, ci = TRUE, ...) : Solution may be nonunique Warning in rq.fit.br(x, y, tau = tau, ci = TRUE, ...) : Solution may be nonunique Warning in rq.fit.br(x, y, tau = tau, ci = TRUE, ...) : Solution may be nonunique Warning in rq.fit.br(x, y, tau = tau, ci = TRUE, ...) : Solution may be nonunique > #Example to illustrate inequality constrained fitting > n <- 100 > p <- 5 > X <- matrix(rnorm(n*p),n,p) > y <- .95*apply(X,1,sum)+rnorm(n) > #constrain slope coefficients to lie between zero and one > R <- cbind(0,rbind(diag(p),-diag(p))) > r <- c(rep(0,p),-rep(1,p)) > rq(y~X,R=R,r=r,method="fnc") Call: rq(formula = y ~ X, method = "fnc", R = R, r = r) Coefficients: (Intercept) X1 X2 X3 X4 X5 -0.3020107 0.7193971 0.9142408 0.9616210 1.0000000 1.0000000 Degrees of freedom: 100 total; 94 residual > > > > cleanEx(); ..nameEx <- "rq.fit.br" > > ### * rq.fit.br > > flush(stderr()); flush(stdout()) > > ### Name: rq.fit.br > ### Title: Quantile Regression Fitting by Exterior Point Methods > ### Aliases: rq.fit.br > ### Keywords: regression > > ### ** Examples > > data(stackloss) > rq.fit.br(stack.x, stack.loss, tau=.73 ,interp=FALSE) $coefficients Air.Flow Water.Temp Acid.Conc. 0.9408138 0.8473490 -0.6308261 $x Air.Flow Water.Temp Acid.Conc. 1 80 27 89 2 80 27 88 3 75 25 90 4 62 24 87 5 62 22 87 6 62 23 87 7 62 24 93 8 62 24 93 9 58 23 87 10 58 18 80 11 58 18 89 12 58 17 88 13 58 18 82 14 58 19 93 15 50 18 89 16 50 18 86 17 50 19 72 18 50 19 79 19 50 20 80 20 56 20 82 21 70 20 91 $y [1] 42 37 37 28 18 18 19 20 15 14 14 13 11 12 8 7 8 8 9 15 15 $residuals [,1] 1 -7.105427e-15 2 -5.630826e+00 3 2.029593e+00 4 4.215043e+00 5 -4.090259e+00 6 -4.937608e+00 7 -1.000000e+00 8 -1.065814e-14 9 -4.174353e+00 10 -5.353391e+00 11 3.240444e-01 12 -4.594328e-01 13 -7.091739e+00 14 -1.243450e-14 15 1.850555e+00 16 -1.041924e+00 17 -9.720838e+00 18 -5.305055e+00 19 -4.521578e+00 20 -2.904809e+00 21 -1.039877e+01 > > > > cleanEx(); ..nameEx <- "rq.fit.sfn" > > ### * rq.fit.sfn > > flush(stderr()); flush(stdout()) > > ### Name: rq.fit.sfn > ### Title: Sparse Regression Quantile Fitting > ### Aliases: rq.fit.sfn > ### Keywords: regression > > ### ** Examples > > n=200 > p=50 > X=rnorm(n*p) > X[abs(X)<2.0]=0 > X=cbind(1,matrix(X,n,p)) > y=.5*apply(X,1,sum)+rnorm(n) > rq.o=rq.fit.sfn(as.matrix.csr(X),y,tmpmax=floor(1e5+exp(-12.1)*(as.matrix.csr(X)@ia[p+1]-1)^2.35)) Loading required package: SparseM > > > > cleanEx(); ..nameEx <- "rq.fit.sfnc" > > ### * rq.fit.sfnc > > flush(stderr()); flush(stdout()) > > ### Name: rq.fit.sfnc > ### Title: Sparse Constrained Regression Quantile Fitting > ### Aliases: rq.fit.sfnc > ### Keywords: regression > > ### ** Examples > > n=200 > p=50 > X=rnorm(n*p) > X[abs(X)<2.0]=0 > X=cbind(1,matrix(X,n,p)) > y=.5*apply(X,1,sum)+rnorm(n) > R=rbind(diag(p+1),-diag(p+1)) > r=c(rep(0,p+1),-rep(1,p+1)) > rq.o=rq.fit.sfnc(as.matrix.csr(X),y,as.matrix.csr(R),r,tmpmax=floor(1e5+exp(-12.1)*(as.matrix.csr(X)@ia[p+1]-1)^2.35)) Loading required package: SparseM > > > > cleanEx(); ..nameEx <- "rqProcess" > > ### * rqProcess > > flush(stderr()); flush(stdout()) > > ### Name: rqProcess > ### Title: Compute Quantile Regression Process > ### Aliases: rqProcess > ### Keywords: regression htest > > ### ** Examples > > > > cleanEx(); ..nameEx <- "rqss" > > ### * rqss > > flush(stderr()); flush(stdout()) > > ### Name: rqss > ### Title: Additive Quantile Regression Smoothing > ### Aliases: rqss rqss.fit [.terms untangle.specials > ### Keywords: regression smooth robust > > ### ** 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) > fit <- rqss(y~qss(x,constraint="N")+z) Loading required package: SparseM > lines(x[-1],fit$coef[1]+fit$coef[-(1:2)]) > fit <- rqss(y~qss(x,constraint="I")+z) > lines(x[-1],fit$coef[1]+fit$coef[-(1:2)],col="blue") > fit <- rqss(y~qss(x,constraint="CI")+z) > lines(x[-1],fit$coef[1]+fit$coef[-(1:2)],col="red") > #Cleanup > rm(list=ls()) > #A bivariate example > data(CobarOre) > fit <- rqss(z~qss(cbind(x,y),lambda=.08),data=CobarOre) Loading required package: tripack > plot(fit) Loading required package: akima > > > > cleanEx(); ..nameEx <- "rrs.test" > > ### * rrs.test > > flush(stderr()); flush(stdout()) > > ### Name: rrs.test > ### Title: Quantile Regression Rankscore Test > ### Aliases: rrs.test > ### Keywords: regression > > ### ** Examples > > # Test that covariates 2 and 3 belong in stackloss model using Wilcoxon scores. > data(stackloss) > rrs.test(stack.x[,1], stack.x[,2:3], stack.loss) $sn [,1] [1,] 11.08975 $ranks [1] 0.23515982 0.02564744 0.38154714 0.47260274 -0.22549020 -0.14862159 [7] -0.08904110 0.13926941 0.02311436 -0.28067541 -0.19696970 -0.31735160 [13] -0.39041096 -0.35388128 0.08967437 -0.15425917 0.20316302 0.27622189 [19] 0.41986749 0.34477124 -0.45433790 > > > > cleanEx(); ..nameEx <- "summary.rq" > > ### * summary.rq > > flush(stderr()); flush(stdout()) > > ### Name: summary.rq > ### Title: Summary methods for Quantile Regression > ### Aliases: summary.rq summary.rqs > ### Keywords: regression > > ### ** Examples > > data(stackloss) > y <- stack.loss > x <- stack.x > summary(rq(y ~ x, method="fn")) # Compute se's for fit using "nid" method. Warning in rq.fit.br(x, y, tau = tau, ci = TRUE, ...) : Solution may be nonunique Call: rq(formula = y ~ x, method = "fn") tau: [1] 0.5 Coefficients: coefficients lower bd upper bd (Intercept) -39.68986 -41.61973 -29.67754 xAir.Flow 0.83188 0.52307 1.14117 xWater.Temp 0.57391 0.32182 1.41090 xAcid.Conc. -0.06087 -0.21348 -0.02891 > summary(rq(y ~ x, ci=FALSE),se="ker") Call: rq(formula = y ~ x, ci = FALSE) tau: [1] 0.5 Coefficients: Value Std. Error t value Pr(>|t|) (Intercept) -39.68986 14.05974 -2.82294 0.01172 xAir.Flow 0.83188 0.24350 3.41632 0.00329 xWater.Temp 0.57391 0.57894 0.99131 0.33543 xAcid.Conc. -0.06087 0.18142 -0.33551 0.74134 > # default "br" alg, and compute kernel method se's > > > > ### *