R Under development (unstable) (2024-05-25 r86622 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

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 an HTML browser interface to help.
Type 'q()' to quit R.

> library(quantreg)
Loading required package: SparseM

Attaching package: 'SparseM'

The following object is masked from 'package:base':

    backsolve

> 
> (dDIR <- system.file("demo", package = "quantreg"))
[1] "D:/RCompile/CRANincoming/R-devel/lib/quantreg/demo"
> set.seed(1) # since some demos randomly generate
> 
> cat("Running demos from package 'quantreg' : \n\n")
Running demos from package 'quantreg' : 

> for(ff in list.files(dDIR, pattern="\\.R$", full.names = TRUE)) {
+    f <- basename(ff)
+    cat("\n", f," :\n", paste(rep.int("-", nchar(f)), collapse=''),
+        "\n", sep='')
+ 
+    source(ff, echo = TRUE)
+ }

Frank.R :
-------

> vFrank <- function(x, df, delta, u) -log(1 - (1 - 
+     exp(-delta))/(1 + exp(-delta * pt(x, df)) * ((1/u) - 1)))/delta

> FrankModel <- function(x, delta, mu, sigma, df, tau) {
+     z <- qt(vFrank(x, df, delta, u = tau), df)
+     mu + sigma * z
+ }

> n <- 200

> df <- 8

> delta <- 8

> set.seed(1989)

> x <- sort(rt(n, df))

> v <- vFrank(x, df, delta, u = runif(n))

> y <- qt(v, df)

> plot(x, y, pch = "o", col = "blue", cex = 0.25)

> Dat <- data.frame(x = x, y = y)

> us <- c(0.25, 0.5, 0.75)

> for (i in 1:length(us)) {
+     v <- vFrank(x, df, delta, u = us[i])
+     lines(x, qt(v, df))
+ }

> cfMat <- matrix(0, 3, length(us))

> trace <- TRUE

> trace <- FALSE

> for (i in 1:length(us)) {
+     tau <- us[i]
+     cat("tau = ", format(tau), ".. ")
+     fit <- nlrq(y ~ FrankModel(x, delta, mu, sigma, df = 8, t .... [TRUNCATED] 
tau =  0.25 .. 
tau =  0.5 .. 
tau =  0.75 .. 

> colnames(cfMat) <- names(coef(fit))

> cfMat
        delta          mu     sigma
[1,] 14.86275 -0.20513700 0.9135934
[2,] 16.65111  0.03327799 0.9578400
[3,] 11.99822  0.11557829 0.9464470

KMvCRQ.R :
--------

> if (requireNamespace("survival", quietly = TRUE)) {
+     n <- 100
+     y <- rchisq(n, 3)
+     c <- rchisq(n, 5)
+     Y <- pmin(y, c)
+     d <-  .... [TRUNCATED] 

MCV.R :
-----

> MCV <- function(lambdas, formula, data, tau = 0.5, 
+     k = 10) {
+     F <- Munge(formula, lambdas = lambdas)
+     f <- rqss(F, data, tau = tau) .... [TRUNCATED] 

> set.seed(1729)

> n <- 200

> x <- sort(runif(n, 0, 20))

> g0 <- function(x, tau) log(x) + 0.2 * (log(x))^3 + 
+     log(x) * qnorm(tau)/4

> y <- g0(x, runif(n))

> D <- data.frame(y = y, x = x)

> lams <- mcvs <- seq(0.02, 5, by = 0.2)

> for (i in 1:length(mcvs)) mcvs[i] <- MCV(lams[i], 
+     y ~ qss(x, lambda = lambdas[1]), D)

> par(mfrow = c(1, 2))

> plot(lams, mcvs, cex = 0.5, lwd = 2, type = "l", xlab = expression(lambda), 
+     ylab = expression(MCV(lambda)))

> lambdastar <- lams[which.min(mcvs)]

> plot(x, y, cex = 0.5, col = "grey")

> f <- rqss(y ~ qss(x, lambda = lambdastar), data = D)

> plot(f, add = TRUE, lwd = 2)

> lines(x, g0(x, 0.5), col = "red", lwd = 2)

> text(10, 1, bquote(lambda == ~.(lambdastar)))

Mammals.R :
---------

> require(quantreg)

> data(Mammals)

> attach(Mammals)

> x <- log(weight)

> xx <- unique(x[order(x)])

> 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 = 0.75)

> taus <- c(0.5, 0.9)

> for (i in 1:length(taus)) {
+     fit <- rqss(y ~ qss(x, lambda = 1, constraint = "C"), tau = taus[i])
+     plot(fit, title = "Running Speed of Mam ..." ... [TRUNCATED] 

> legend(4, 2, c("Median", "0.9 Quantile"), lty = 1, 
+     col = 1:2, lwd = 1.5)

> plot(fit, title = "Running Speed of Mammals", band = "both", 
+     col = i, lwd = 1.5)

> xy <- fit$qss[[1]]$xyz

> xx <- xy[, 1]

> yhat <- fit$coef[1] + xy[, 2]

> g <- diff(yhat)/diff(xx)

> plot(xx[-1], g, main = "Fitted Slopes of Running Speed", 
+     xlab = "Weight in log(Kg)", ylab = "dlog(Speed) /dlog(Weight)")

Mel.R :
-----

> require(splines)
Loading required package: splines

> if (interactive()) {
+     oldpar <- par(ask = TRUE)
+     data(MelTemp)
+     x <- MelTemp[-3650]
+     y <- MelTemp[-1]
+     s <- (x < 40)
+      .... [TRUNCATED] 

Mel2.R :
------

> require(splines)

> if (interactive()) {
+     oldpar <- par(ask = TRUE)
+     data(MelTemp)
+     x <- MelTemp[-3650]
+     y <- MelTemp[-1]
+     s <- (x < 40)
+      .... [TRUNCATED] 

Panel.R :
-------

> library(quantreg)

> rq.fit.panel <- function(X, y, s, w = c(1/3, 1/3, 
+     1/3), taus = c(0.25, 0.5, 0.75), lambda = 0) {
+     require(SparseM)
+     require(quantre .... [TRUNCATED] 

> n <- 3

> T <- 50

> nT <- n * T

> u1 <- rnorm(T)

> u2 <- rnorm(T)

> u3 <- rnorm(T)

> x1 <- rnorm(T, 1, 0.85)

> x2 <- rnorm(T, 4, 1)

> x3 <- rnorm(T, 7, 1)

> beta1 <- 1

> beta2 <- 0.5

> y1 <- 0 + beta1 * x1 + (beta2 * x1) * u1

> y2 <- 4 + beta1 * x2 + (beta2 * x2) * u2

> y3 <- 8 + beta1 * x3 + (beta2 * x3) * u3

> plot(c(0, 9), c(0, 25), type = "n", xlab = expression(x[it]), 
+     ylab = expression(y[it]))

> points(x1, y1, pch = 15)

> points(x2, y2, pch = 15, col = "blue")

> points(x3, y3, pch = 15, col = "red")

> legend(1, 17, paste("i = ", 1:3, sep = ""), pch = 15, 
+     col = c("black", "blue", "red"))

> ya <- c(y1, y2, y3)

> xa <- c(x1, x2, x3)

> taus <- c(0.25, 0.5, 0.75)

> xx <- seq(min(xa), max(xa), 0.25)

> f <- coef(rq(ya ~ xa, tau = taus))

> yy <- cbind(1, xx) %*% f

> for (i in 1:3) lines(xx, yy[, i], col = "grey")

> s <- rep(1:n, rep(T, n))

> fp <- rq.fit.panel(xa, ya, s)$coef

> bhat <- fp[1:3]

> fehat <- fp[4:6]

> xx1 <- seq(min(x1), max(x1), 0.25)

> for (i in 1:3) {
+     yy1 <- fehat[1] + bhat[i] * xx1
+     lines(xx1, yy1, col = "black")
+ }

> xx2 <- seq(min(x2), max(x2), 0.25)

> for (i in 1:3) {
+     yy2 <- fehat[2] + bhat[i] * xx2
+     lines(xx2, yy2, col = "blue")
+ }

> xx3 <- seq(min(x3), max(x3), 0.25)

> for (i in 1:3) {
+     yy3 <- fehat[3] + bhat[i] * xx3
+     lines(xx3, yy3, col = "red")
+ }

Polson.R :
--------

> dgp <- function(n) {
+     x <- 0:n/n
+     y <- rnorm(n + 1, 5 * sin(2 * pi * x), 0.5 + exp(1.5 * sin(4 * 
+         pi * x)))
+     data.frame(x = .... [TRUNCATED] 

> D <- dgp(1000)

> plot(D$x, D$y, cex = 0.5)

> taus <- 1:9/10

> for (i in 1:length(taus)) plot(rqss(y ~ qss(x, lambda = 1/10), 
+     tau = taus[i], data = D), rug = FALSE, add = TRUE)

RB-r.R :
------

> if (require(MASS)) {
+     require(quantreg)
+     X <- cbind(1, matrix(rnorm(500), 100, 5))
+     y <- rnorm(100)
+     R <- matrix(rnorm(18), 3, 6 .... [TRUNCATED] 
Loading required package: MASS

arqss.R :
-------

> n <- 2000

> x <- 1:n/n

> noise <- rgamma(n, 3, 1)

> g0 <- function(x) sin(10 * x)

> y <- g0(x) + noise

> arqss <- function(x, y, tau, g0 = NULL) {
+     g <- function(lam, y, x, tau) AIC(rqss(y ~ qss(x, lambda = lam), 
+         tau = tau), k = -1)
+    .... [TRUNCATED] 

> arqss(x, y, 0.5, g0)

cobar.R :
-------

> require(quantreg)

> if (requireNamespace("interp")) {
+     do.rgl <- interactive() && require(rgl)
+     data(CobarOre)
+     fit <- rqss(z ~ qss(cbind(x, y), lambda = .... [TRUNCATED] 
Loading required namespace: interp

combos.R :
--------

> H <- combos(20, 3)

> if (!require("rgl", quietly = TRUE)) {
+     warning("The package rgl is needed for plotting")
+ } else {
+     if (interactive()) {
+         plot3 .... [TRUNCATED] 

cpoint.R :
--------

> x <- runif(200, 0, 10)

> u <- rnorm(200)/5

> y <- 1 + x - 0.5 * (x - 3) * (x > 3) + u

> plot(y ~ x, cex = 0.5, col = "grey")

> z <- rqss(y ~ qss(x, lambda = 10), tau = 0.5)

> plot(z, col = "dark blue")

> eps <- 1e-05

> Nz <- abs(z$resid[1:200]) < eps

> Nj <- abs(z$resid[201:398]) > eps

> xx <- z$qss[[1]]$xyz[, 1]

> yy <- z$coef[1] + z$qss[[1]]$xyz[, 2]

> xj <- xx[3:200]

> yj <- yy[3:200]

> points(xx[Nz], yy[Nz], col = "green")

> points(xj[Nj], yj[Nj], col = "red")

> print(paste("Number of zero residuals:  ", sum(Nz)))
[1] "Number of zero residuals:   4"

> print(paste("Number of jumps in slope:  ", sum(Nj)))
[1] "Number of jumps in slope:   2"

> legend(6, 3, c("Derivative Jumps", "Zero residuals"), 
+     pch = "o", col = c("red", "green"))

crquis.R :
--------

> if (requireNamespace("survival", quietly = TRUE)) {
+     data(uis)
+     Surv <- survival::Surv
+     fit <- crq(Surv(log(TIME), CENSOR) ~ ND1 + ND .... [TRUNCATED] 

engel1.R :
--------

> data(engel)

> plot(foodexp ~ income, data = engel, cex = 0.5, col = "blue", 
+     xlab = "Household Income", ylab = "Food Expenditure")

> z <- rq(foodexp ~ income, tau = 0.5, data = engel)

> abline(z, col = "dark blue")

> abline(lm(foodexp ~ income, data = engel), lty = 2, 
+     col = "red")

> taus <- c(0.05, 0.1, 0.25, 0.75, 0.9, 0.95)

> nt <- length(taus)

> for (i in 1:length(taus)) {
+     abline(rq(foodexp ~ income, tau = taus[i], data = engel), 
+         col = "gray")
+ }

> legend("bottomright", c("L1 (tau = .50)", "OLS", paste("tau= ", 
+     formatC(rev(taus)))), col = c("dark blue", "red", rep("gray", 
+     nt)), lt .... [TRUNCATED] 

engel2.R :
--------

> data(engel)

> x.poor <- quantile(engel[, "income"], 0.1)

> x.rich <- quantile(engel[, "income"], 0.9)

> z <- rq(foodexp ~ income, tau = -1, data = engel)

> ps <- z$sol["tau", ]

> coefs <- z$sol[4:5, ]

> qs.poor <- c(c(1, x.poor) %*% coefs)

> qs.rich <- c(c(1, x.rich) %*% coefs)

> par(mfrow = c(1, 2))

> plot(c(ps, ps), c(qs.poor, qs.rich), type = "n", xlab = expression(tau), 
+     ylab = "quantile")

> plot(stepfun(ps, c(qs.poor[1], qs.poor)), do.points = FALSE, 
+     add = TRUE)

> plot(stepfun(ps, c(qs.poor[1], qs.rich)), do.points = FALSE, 
+     add = TRUE, col.hor = "gray", col.vert = "gray")

> ps.wts <- (c(0, diff(ps)) + c(diff(ps), 0))/2

> ap <- akj(qs.poor, z = qs.poor, p = ps.wts)

> ar <- akj(qs.rich, z = qs.rich, p = ps.wts)

> plot(c(qs.poor, qs.rich), c(ap$dens, ar$dens), type = "n", 
+     xlab = "Food Expenditure", ylab = "Density")

> lines(qs.rich, ar$dens, col = "gray")

> lines(qs.poor, ap$dens, col = "black")

> legend("topright", c("poor", "rich"), lty = c(1, 1), 
+     col = c("black", "gray"))

hinged.R :
--------

> require(quantreg)

> if (requireNamespace("interp")) {
+     do.rgl <- interactive() && require(rgl)
+     n <- 1000
+     x <- runif(n)
+     y <- runif(n)
+     z <- - .... [TRUNCATED] 
NULL

panelfig.R :
----------

> library(quantreg)

> rq.fit.panel <- function(X, y, s, w = c(1/3, 1/3, 
+     1/3), taus = c(0.25, 0.5, 0.75), lambda = 0) {
+     require(SparseM)
+     require(quantre .... [TRUNCATED] 

> n <- 3

> T <- 50

> nT <- n * T

> u1 <- rnorm(T)

> u2 <- rnorm(T)

> u3 <- rnorm(T)

> x1 <- rnorm(T, 1, 0.85)

> x2 <- rnorm(T, 4, 1)

> x3 <- rnorm(T, 7, 1)

> beta1 <- 1

> beta2 <- 0.5

> y1 <- 0 + beta1 * x1 + (beta2 * x1) * u1

> y2 <- 4 + beta1 * x2 + (beta2 * x2) * u2

> y3 <- 8 + beta1 * x3 + (beta2 * x3) * u3

> plot(c(0, 9), c(0, 25), type = "n", xlab = expression(x[it]), 
+     ylab = expression(y[it]))

> points(x1, y1, pch = 15)

> points(x2, y2, pch = 15, col = "blue")

> points(x3, y3, pch = 15, col = "red")

> legend(1, 17, paste("i = ", 1:3, sep = ""), pch = 15, 
+     col = c("black", "blue", "red"))

> ya <- c(y1, y2, y3)

> xa <- c(x1, x2, x3)

> taus <- c(0.25, 0.5, 0.75)

> xx <- seq(min(xa), max(xa), 0.25)

> f <- coef(rq(ya ~ xa, tau = taus))

> yy <- cbind(1, xx) %*% f

> for (i in 1:3) lines(xx, yy[, i], col = "grey")

> s <- rep(1:n, rep(T, n))

> fp <- rq.fit.panel(xa, ya, s)$coef

> bhat <- fp[1:3]

> fehat <- fp[4:6]

> xx1 <- seq(min(x1), max(x1), 0.25)

> for (i in 1:3) {
+     yy1 <- fehat[1] + bhat[i] * xx1
+     lines(xx1, yy1, col = "black")
+ }

> xx2 <- seq(min(x2), max(x2), 0.25)

> for (i in 1:3) {
+     yy2 <- fehat[2] + bhat[i] * xx2
+     lines(xx2, yy2, col = "blue")
+ }

> xx3 <- seq(min(x3), max(x3), 0.25)

> for (i in 1:3) {
+     yy3 <- fehat[3] + bhat[i] * xx3
+     lines(xx3, yy3, col = "red")
+ }

predemo.R :
---------

> par(ask = TRUE)

> sub <- "Classical Gaussian Version"

> plot(cars$speed, cars$dist, xlim = c(0, 27), ylim = c(-10, 
+     130), main = "Ezekiel Data", sub = sub, xlab = "speed (mph)", 
+     ylab = "stopp ..." ... [TRUNCATED] 

> ss = seq(0, 27, len = 100)

> flm <- lm(dist ~ speed + I(speed^2), data = cars)

> cflm <- predict(flm, list(speed = ss), data = cars, 
+     interval = "confidence")

> pflm <- predict(flm, list(speed = ss), data = cars, 
+     interval = "prediction")

> matlines(ss, pflm, lty = c(1, 2, 2), col = c("black", 
+     "red", "red"))

> matlines(ss, cflm, lty = c(1, 2, 2), col = c("black", 
+     "blue", "blue"))

> fl = rq(dist ~ speed + I(speed^2), data = cars, tau = 0.25)

> fu = rq(dist ~ speed + I(speed^2), data = cars, tau = 0.75)

> sub <- "Covariance matrix estimation using nid"

> plot(cars$speed, cars$dist, xlim = c(0, 27), ylim = c(-10, 
+     130), main = "Ezekiel Data", sub = sub, xlab = "speed (mph)", 
+     ylab = "stopp ..." ... [TRUNCATED] 

> pl <- predict(fl, list(speed = ss), interval = "confidence", 
+     se = "nid")

> pu <- predict(fu, list(speed = ss), interval = "confidence", 
+     se = "nid")

> matlines(ss, pl, lty = c(1, 2, 0), col = c("black", 
+     "red", "red"))

> matlines(ss, pu, lty = c(1, 0, 2), col = c("black", 
+     "red", "red"))

> sub <- "Covariance matrix estimation using xy-bootstrap"

> plot(cars$speed, cars$dist, xlim = c(0, 27), ylim = c(-10, 
+     130), main = "Ezekiel Data", sub = sub, xlab = "speed (mph)", 
+     ylab = "stopp ..." ... [TRUNCATED] 

> pl <- predict(fl, list(speed = ss), interval = "confidence", 
+     se = "boot", bsmethod = "xy")

> pu <- predict(fu, list(speed = ss), interval = "confidence", 
+     se = "boot", bsmethod = "xy")

> matlines(ss, pl, lty = c(1, 2, 0), col = c("black", 
+     "red", "red"))

> matlines(ss, pu, lty = c(1, 0, 2), col = c("black", 
+     "red", "red"))

> sub <- "Covariance matrix estimation using weighted xy-bootstrap"

> plot(cars$speed, cars$dist, xlim = c(0, 27), ylim = c(-10, 
+     130), main = "Ezekiel Data", sub = sub, xlab = "speed (mph)", 
+     ylab = "stopp ..." ... [TRUNCATED] 

> pl <- predict(fl, list(speed = ss), interval = "confidence", 
+     se = "boot", bsmethod = "wxy")

> pu <- predict(fu, list(speed = ss), interval = "confidence", 
+     se = "boot", bsmethod = "wxy")

> matlines(ss, pl, lty = c(1, 2, 0), col = c("black", 
+     "red", "red"))

> matlines(ss, pu, lty = c(1, 0, 2), col = c("black", 
+     "red", "red"))

> sub <- "Percentile method using xy-bootstrap"

> plot(cars$speed, cars$dist, xlim = c(0, 27), ylim = c(-10, 
+     130), main = "Ezekiel Data", sub = sub, xlab = "speed (mph)", 
+     ylab = "stopp ..." ... [TRUNCATED] 

> pl <- predict(fl, list(speed = ss), interval = "confidence", 
+     type = "percentile", se = "boot", bsmethod = "xy")

> pu <- predict(fu, list(speed = ss), interval = "confidence", 
+     type = "percentile", se = "boot", bsmethod = "xy")

> matlines(ss, pl, lty = c(1, 2, 0), col = c("black", 
+     "red", "red"))

> matlines(ss, pu, lty = c(1, 0, 2), col = c("black", 
+     "red", "red"))

> sub <- "Direct method using iid covariance matrix"

> plot(cars$speed, cars$dist, xlim = c(0, 27), ylim = c(-10, 
+     130), main = "Ezekiel Data", sub = sub, xlab = "speed (mph)", 
+     ylab = "stopp ..." ... [TRUNCATED] 

> pl <- predict(fl, list(speed = ss), interval = "confidence", 
+     type = "direct", se = "iid")

> pu <- predict(fu, list(speed = ss), interval = "confidence", 
+     type = "direct", se = "iid")

> matlines(ss, pl, lty = c(1, 2, 0), col = c("black", 
+     "red", "red"))

> matlines(ss, pu, lty = c(1, 0, 2), col = c("black", 
+     "red", "red"))

> sub <- "Direct method using nid covariance matrix"

> plot(cars$speed, cars$dist, xlim = c(0, 27), ylim = c(-10, 
+     130), main = "Ezekiel Data", sub = sub, xlab = "speed (mph)", 
+     ylab = "stopp ..." ... [TRUNCATED] 

> pl <- predict(fl, list(speed = ss), interval = "confidence", 
+     type = "direct", se = "nid")

> pu <- predict(fu, list(speed = ss), interval = "confidence", 
+     type = "direct", se = "nid")

> matlines(ss, pl, lty = c(1, 2, 0), col = c("black", 
+     "red", "red"))

> matlines(ss, pu, lty = c(1, 0, 2), col = c("black", 
+     "red", "red"))

> par(ask = FALSE)

rqsslasso.R :
-----------

> n <- 100

> p <- 9

> q <- 3

> beta <- c(rep(1, q), rep(0, p - q))

> w <- matrix(rnorm(n * p), n, p)

> x <- runif(n, 0, 10)

> z <- runif(n, 0, 10)

> y <- w %*% beta + sin(x) + (z^2)/50 + rnorm(n)/5

> d <- data.frame(w, x, y, z)

> f <- rqss(y ~ w + qss(x, lambda = 3) + qss(z, lambda = 2), 
+     method = "lasso", lambda = 3, data = d)

> plot(f, bands = "both", bcol = c("lightsteelblue", 
+     "lightsteelblue4"))

stack.R :
-------

> require(quantreg)

> data(stackloss)

> logLik.rq.process <- function(fit) {
+     y <- model.response(model.frame(fit))
+     fhat <- predict(fit, type = "fhat")
+     fy <- mapply(functi .... [TRUNCATED] 

> f0 <- rq(stack.loss ~ 1, tau = -1)

> f1 <- rq(stack.loss ~ stack.x, tau = -1)

> l0 <- logLik(f0)

> l1 <- logLik(f1)

> f0 <- rq(stack.loss ~ 1, tau = 1:19/20)

> f1 <- rq(stack.loss ~ stack.x, tau = 1:19/20)

> l0 <- logLik(f0)

> l1 <- logLik(f1)

subset.R :
--------

> require(quantreg)

> n <- 200

> x <- sort(rchisq(n, 4))

> z <- rnorm(n)

> s <- sample(1:n, n/2)

> y <- log(x) + rnorm(n)/5

> D = data.frame(y = y, x = x, z = z, s = (1:n) %in% 
+     s)

> plot(x, y, cex = 0.5, col = "grey")

> points(x[s], y[s], col = "pink", cex = 0.5)

> lam = 0.2

> f0 <- rqss(y ~ qss(x, lambda = lam) + z, subset = s)

> f1 <- rqss(y ~ qss(x, lambda = lam) + z, subset = s, 
+     data = D)

> plot(f0, add = TRUE, col = 2, lwd = 3)

> plot(f1, add = TRUE, col = 4, lwd = 3)
Warning messages:
1: In summary.rq(x, se = se, R = R, covariance = TRUE) :
  5 non-positive fis
2: In survfit.coxph(x) :
  the model contains interactions; the default curve based on columm means of the X matrix is almost certainly not useful. Consider adding a newdata argument.
3: In pt(abs(ptab[, 3]), n - edf) : NaNs produced
4: In pf(ntab[i - 1, 4], ntab[i - 1, 1], n - edf) : NaNs produced
5: In summary.rq(object, cov = TRUE, ...) : 2 non-positive fis
6: In summary.rq(object, cov = TRUE, ...) : 2 non-positive fis
> 
> 
> proc.time()
   user  system elapsed 
  33.21    5.68   39.29