piqr: Penalized Quantile Regression Coefficients Modeling
In qrcmNP: Nonlinear and Penalized Quantile Regression Coefficients Modeling
piqr R Documentation
Penalized Quantile Regression Coefficients Modeling
Description
This package implements a penalized Frumento and Bottai's (2016) method for quantile regression coefficient modeling (qrcm), in which quantile regression coefficients are described by (flexible) parametric functions of the order of the quantile.
This package fits lasso qrcm using pathwise coordinate descent algorithm.
Usage
piqr(formula, formula.p = ~ slp(p, 3), weights, data, s, nlambda=100,
lambda.min.ratio=ifelse(nobs<nvars, 0.01, 0.0001), lambda,
tol=1e-6, maxit=100, display=TRUE)
Arguments
formula
a two-sided formula of the form y ~ x1 + x2 + ...:
a symbolic description of the quantile regression model.
formula.p
a one-sided formula of the form ~ b1(p, ...) + b2(p, ...) + ..., describing how
quantile regression coefficients depend on p, the order of the quantile.
weights
an optional vector of weights to be used in the fitting process.
data
an optional data frame, list or environment containing the variables in formula.
s
an optional 0/1 matrix that permits excluding some model coefficients
(see ‘Examples’).
nlambda
the number of lambda values - default is 100.
lambda.min.ratio
Smallest value for lambda, as a fraction of lambda.max. The default depends on the sample size nobs relative to the number of variables nvars. If nobs > nvars, the default is 0.0001, close to zero. If nobs < nvars, the default is 0.01.
lambda
A user supplied lambda sequence.
display
if TRUE something is printed - default is TRUE.
tol
convergence criterion for numerical optimization - default is 1e-6.
maxit
maximum number of iterations - default is 100.
Details
Quantile regression permits modeling conditional quantiles of a response variabile,
given a set of covariates. A linear model is used to describe the conditional
quantile function:
Q(p | x) = \beta_0(p) + \beta_1(p)x_1 + \beta_2(p)x_2 + \ldots.
The model coefficients \beta(p) describe the effect of covariates on the p-th
quantile of the response variable. Usually, one or more
quantiles are estimated, corresponding to different values of p.
Assume that each coefficient can be expressed as a parametric function of p of the form:
\beta(p | \theta) = \theta_{0} + \theta_1 b_1(p) + \theta_2 b_2(p) + \ldots
where b_1(p), b_2(p, \ldots) are known functions of p.
If q is the dimension of
x = (1, x_1, x_2, \ldots)
and k is that of
b(p) = (1, b_1(p), b_2(p), \ldots),
the entire conditional quantile function is described by a
q \times k matrix \theta of model parameters.
Users are required to specify two formulas: formula describes the regression model,
while formula.p identifies the 'basis' b(p).
By default, formula.p = ~ slp(p, k = 3), a 3rd-degree shifted
Legendre polynomial (see slp). Any user-defined function b(p, \ldots)
can be used, see ‘Examples’.
Estimation of penalized \theta is carried out by minimizing a penalized integrated loss function,
corresponding to the integral, over p, of the penalized loss function of standard quantile regression. This
motivates the acronym piqr (penalized integrated quantile regression).
See details in iqr
Value
An object of class “piqr”, a list containing the following items:
call
the matched call.
lambda
The actual sequence of lambda values used.
coefficients
a list of estimated model parameters describing the fitted quantile function along the path.
minimum
the value of the minimized integrated loss function for each value of lambda.
dl
a matrix of gradient values along the path.
df
The number of nonzero coefficients for each value of lambda.
seqS
a list containg each matrix s for each value of lambda.
internal
a list containing some initial object.
Note
By expressing quantile regression coefficients as functions of p, a parametric model for the conditional
quantile function is specified. The induced PDF and CDF can be used as diagnostic tools.
Negative values of PDF indicate quantile crossing, i.e., the conditional quantile function is not
monotonically increasing. Null values of PDF indicate observations that lie outside the
estimated support of the data, defined by quantiles of order 0 and 1. If null or negative PDF
values occur for a relatively large proportion of data, the model is probably misspecified or ill-defined.
If the model is correct, the fitted CDF should approximately follow a Uniform(0,1) distribution.
This idea is used to implement a goodness-of-fit test, see summary.iqr
and test.fit.
The intercept can be excluded from formula, e.g.,
iqr(y ~ -1 + x). This, however, implies that when x = 0,
y is always 0. See example 5 in ‘Examples’.
The intercept can also be removed from formula.p.
This is recommended if the data are bounded. For example, for strictly positive data,
use iqr(y ~ 1, formula.p = -1 + slp(p,3)) to force the smallest quantile
to be zero. See example 6 in ‘Examples’.
Author(s)
Gianluca Sottile gianluca.sottile@unipa.it
References
Sottile G, Frumento P, Chiodi M, Bottai M. (2020). A penalized approach to covariate selection through quantile regression coefficient models. Statistical Modelling, 20(4), pp 369-385. doi:10.1177/1471082X19825523.
Frumento, P., and Bottai, M. (2016). Parametric modeling of quantile regression coefficient functions. Biometrics, 72(1), pp 74-84, doi:10.1111/biom.12410.
Friedman, J., Hastie, T. and Tibshirani, R. (2008). Regularization Paths for Generalized Linear Models via Coordinate Descent. Journal of Statistical Software, Vol. 33(1), pp 1-22 Feb 2010.
See Also
summary.piqr, plot.piqr, predict.piqr,
for summary, plotting, and prediction.
gof.piqr to select the best value of the tuning parameter though AIC, BIC, GIC, GCV criteria.
Examples
##### Using simulated data in all examples
##### Example 1
set.seed(1234)
n <- 300
x1 <- rexp(n)
x2 <- runif(n, 0, 5)
x <- cbind(x1,x2)
b <- function(p){matrix(cbind(1, qnorm(p), slp(p, 2)), nrow=4, byrow=TRUE)}
theta <- matrix(0, nrow=3, ncol=4); theta[, 1] <- 1; theta[1,2] <- 1; theta[2:3,3] <- 2
qy <- function(p, theta, b, x){rowSums(x * t(theta %*% b(p)))}
y <- qy(runif(n), theta, b, cbind(1, x))
s <- matrix(1, nrow=3, ncol=4); s[1,3:4] <- 0; s[2:3, 2] <- 0
obj <- piqr(y ~ x1 + x2, formula.p = ~ I(qnorm(p)) + slp(p, 2), s=s, nlambda=50)
best <- gof.piqr(obj, method="AIC", plot=FALSE)
best2 <- gof.piqr(obj, method="BIC", plot=FALSE)
summary(obj, best$posMinLambda)
summary(obj, best2$posMinLambda)
## Not run:
##### other examples
set.seed(1234)
n <- 1000
q <- 5
k <- 3
X <- matrix(abs(rnorm(n*q)), n, q)
rownames(X) <- 1:n
colnames(X) <- paste0("X", 1:q)
theta <- matrix(c(3, 1.5, 1, 1,
2, 1, 1, 1,
0, 0, 0, 0,
0, 0, 0, 0,
1.5, 1, 1, 1,
0, 0, 0, 0),
ncol=(k+1), byrow=TRUE)
rownames(theta) <- c("(intercept)", paste0("X", 1:q))
colnames(theta) <- c("(intercept)", "slp(p,1)", "slp(p,2)", "slp(p,3)")
B <- function(p, k){matrix(cbind(1, slp(p, k)), nrow=(k+1), byrow=TRUE)}
Q <- function(p, theta, B, k, X){rowSums(X * t(theta %*% B(p, k)))}
pp <- runif(n)
y <- Q(p=pp, theta=theta, B=B, k=k, X=cbind(1, X))
m1 <- piqr(y ~ X, formula.p = ~ slp(p, k))
best1 <- gof.piqr(m1, method="AIC", plot=FALSE)
best2 <- gof.piqr(m1, method="BIC", plot=FALSE)
summary(m1, best1$posMinLambda)
summary(m1, best2$posMinLambda)
par(mfrow = c(1,3)); plot(m1, xvar="lambda");
plot(m1, xvar="objective"); plot(m1, xvar="grad")
set.seed(1234)
n <- 1000
q <- 6
k <- 4
# x <- runif(n)
X <- matrix(abs(rnorm(n*q)), n, q)
rownames(X) <- 1:n
colnames(X) <- paste0("X", 1:q)
theta <- matrix(c(1, 2, 0, 0, 0,
2, 0, 1, 0, 0,
0, 0, 0, 0, 0,
1, 0, 0, 1, -1.2,
0, 0, 0, 0, 0,
1.5, 0, .5, 0, 0,
0, 0, 0, 0, 0),
ncol=(k+1), byrow=TRUE)
rownames(theta) <- c("(intercept)", paste0("X", 1:q))
colnames(theta) <- c("(intercept)", "qnorm(p)", "p", "log(p)", "log(1-p)")
B <- function(p, k){matrix(cbind(1, qnorm(p), p, log(p), log(1-p)), nrow=(k+1), byrow=TRUE)}
Q <- function(p, theta, B, k, X){rowSums(X * t(theta %*% B(p, k)))}
pp <- runif(n)
y <- Q(p=pp, theta=theta, B=B, k=k, X=cbind(1, X))
s <- matrix(1, q+1, k+1); s[2:(q+1), 2] <- 0; s[1, 3:(k+1)] <- 0; s[2:3, 4:5] <- 0
s[4:5, 3] <- 0; s[6:7, 4:5] <- 0
m2 <- piqr(y ~ X, formula.p = ~ qnorm(p) + p + I(log(p)) + I(log(1-p)), s=s)
best1 <- gof.piqr(m2, method="AIC", plot=FALSE)
best2 <- gof.piqr(m2, method="BIC", plot=FALSE)
summary(m2, best1$posMinLambda)
summary(m2, best2$posMinLambda)
par(mfrow = c(1,3)); plot(m2, xvar="lambda");
plot(m2, xvar="objective"); plot(m2, xvar="grad")
## End(Not run)
# see the documentation for 'summary.piqr', and 'plot.piqr'
qrcmNP documentation built on May 29, 2024, 8:29 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
| piqr | R Documentation |
Penalized Quantile Regression Coefficients Modeling
Description
This package implements a penalized Frumento and Bottai's (2016) method for quantile regression coefficient modeling (qrcm), in which quantile regression coefficients are described by (flexible) parametric functions of the order of the quantile. This package fits lasso qrcm using pathwise coordinate descent algorithm.
Usage
piqr(formula, formula.p = ~ slp(p, 3), weights, data, s, nlambda=100,
lambda.min.ratio=ifelse(nobs<nvars, 0.01, 0.0001), lambda,
tol=1e-6, maxit=100, display=TRUE)
Arguments
formula |
a two-sided formula of the form |
formula.p |
a one-sided formula of the form |
weights |
an optional vector of weights to be used in the fitting process. |
data |
an optional data frame, list or environment containing the variables in |
s |
an optional 0/1 matrix that permits excluding some model coefficients (see ‘Examples’). |
nlambda |
the number of lambda values - default is 100. |
lambda.min.ratio |
Smallest value for lambda, as a fraction of lambda.max. The default depends on the sample size nobs relative to the number of variables nvars. If nobs > nvars, the default is 0.0001, close to zero. If nobs < nvars, the default is 0.01. |
lambda |
A user supplied lambda sequence. |
display |
if TRUE something is printed - default is TRUE. |
tol |
convergence criterion for numerical optimization - default is 1e-6. |
maxit |
maximum number of iterations - default is 100. |
Details
Quantile regression permits modeling conditional quantiles of a response variabile, given a set of covariates. A linear model is used to describe the conditional quantile function:
Q(p | x) = \beta_0(p) + \beta_1(p)x_1 + \beta_2(p)x_2 + \ldots.
The model coefficients \beta(p) describe the effect of covariates on the p-th
quantile of the response variable. Usually, one or more
quantiles are estimated, corresponding to different values of p.
Assume that each coefficient can be expressed as a parametric function of p of the form:
\beta(p | \theta) = \theta_{0} + \theta_1 b_1(p) + \theta_2 b_2(p) + \ldots
where b_1(p), b_2(p, \ldots) are known functions of p.
If q is the dimension of
x = (1, x_1, x_2, \ldots)
and k is that of
b(p) = (1, b_1(p), b_2(p), \ldots),
the entire conditional quantile function is described by a
q \times k matrix \theta of model parameters.
Users are required to specify two formulas: formula describes the regression model,
while formula.p identifies the 'basis' b(p).
By default, formula.p = ~ slp(p, k = 3), a 3rd-degree shifted
Legendre polynomial (see slp). Any user-defined function b(p, \ldots)
can be used, see ‘Examples’.
Estimation of penalized \theta is carried out by minimizing a penalized integrated loss function,
corresponding to the integral, over p, of the penalized loss function of standard quantile regression. This
motivates the acronym piqr (penalized integrated quantile regression).
See details in iqr
Value
An object of class “piqr”, a list containing the following items:
call |
the matched call. |
lambda |
The actual sequence of lambda values used. |
coefficients |
a list of estimated model parameters describing the fitted quantile function along the path. |
minimum |
the value of the minimized integrated loss function for each value of lambda. |
dl |
a matrix of gradient values along the path. |
df |
The number of nonzero coefficients for each value of lambda. |
seqS |
a list containg each matrix s for each value of lambda. |
internal |
a list containing some initial object. |
Note
By expressing quantile regression coefficients as functions of p, a parametric model for the conditional
quantile function is specified. The induced PDF and CDF can be used as diagnostic tools.
Negative values of PDF indicate quantile crossing, i.e., the conditional quantile function is not
monotonically increasing. Null values of PDF indicate observations that lie outside the
estimated support of the data, defined by quantiles of order 0 and 1. If null or negative PDF
values occur for a relatively large proportion of data, the model is probably misspecified or ill-defined.
If the model is correct, the fitted CDF should approximately follow a Uniform(0,1) distribution.
This idea is used to implement a goodness-of-fit test, see summary.iqr
and test.fit.
The intercept can be excluded from formula, e.g.,
iqr(y ~ -1 + x). This, however, implies that when x = 0,
y is always 0. See example 5 in ‘Examples’.
The intercept can also be removed from formula.p.
This is recommended if the data are bounded. For example, for strictly positive data,
use iqr(y ~ 1, formula.p = -1 + slp(p,3)) to force the smallest quantile
to be zero. See example 6 in ‘Examples’.
Author(s)
Gianluca Sottile gianluca.sottile@unipa.it
References
Sottile G, Frumento P, Chiodi M, Bottai M. (2020). A penalized approach to covariate selection through quantile regression coefficient models. Statistical Modelling, 20(4), pp 369-385. doi:10.1177/1471082X19825523.
Frumento, P., and Bottai, M. (2016). Parametric modeling of quantile regression coefficient functions. Biometrics, 72(1), pp 74-84, doi:10.1111/biom.12410.
Friedman, J., Hastie, T. and Tibshirani, R. (2008). Regularization Paths for Generalized Linear Models via Coordinate Descent. Journal of Statistical Software, Vol. 33(1), pp 1-22 Feb 2010.
See Also
summary.piqr, plot.piqr, predict.piqr,
for summary, plotting, and prediction.
gof.piqr to select the best value of the tuning parameter though AIC, BIC, GIC, GCV criteria.
Examples
##### Using simulated data in all examples
##### Example 1
set.seed(1234)
n <- 300
x1 <- rexp(n)
x2 <- runif(n, 0, 5)
x <- cbind(x1,x2)
b <- function(p){matrix(cbind(1, qnorm(p), slp(p, 2)), nrow=4, byrow=TRUE)}
theta <- matrix(0, nrow=3, ncol=4); theta[, 1] <- 1; theta[1,2] <- 1; theta[2:3,3] <- 2
qy <- function(p, theta, b, x){rowSums(x * t(theta %*% b(p)))}
y <- qy(runif(n), theta, b, cbind(1, x))
s <- matrix(1, nrow=3, ncol=4); s[1,3:4] <- 0; s[2:3, 2] <- 0
obj <- piqr(y ~ x1 + x2, formula.p = ~ I(qnorm(p)) + slp(p, 2), s=s, nlambda=50)
best <- gof.piqr(obj, method="AIC", plot=FALSE)
best2 <- gof.piqr(obj, method="BIC", plot=FALSE)
summary(obj, best$posMinLambda)
summary(obj, best2$posMinLambda)
## Not run:
##### other examples
set.seed(1234)
n <- 1000
q <- 5
k <- 3
X <- matrix(abs(rnorm(n*q)), n, q)
rownames(X) <- 1:n
colnames(X) <- paste0("X", 1:q)
theta <- matrix(c(3, 1.5, 1, 1,
2, 1, 1, 1,
0, 0, 0, 0,
0, 0, 0, 0,
1.5, 1, 1, 1,
0, 0, 0, 0),
ncol=(k+1), byrow=TRUE)
rownames(theta) <- c("(intercept)", paste0("X", 1:q))
colnames(theta) <- c("(intercept)", "slp(p,1)", "slp(p,2)", "slp(p,3)")
B <- function(p, k){matrix(cbind(1, slp(p, k)), nrow=(k+1), byrow=TRUE)}
Q <- function(p, theta, B, k, X){rowSums(X * t(theta %*% B(p, k)))}
pp <- runif(n)
y <- Q(p=pp, theta=theta, B=B, k=k, X=cbind(1, X))
m1 <- piqr(y ~ X, formula.p = ~ slp(p, k))
best1 <- gof.piqr(m1, method="AIC", plot=FALSE)
best2 <- gof.piqr(m1, method="BIC", plot=FALSE)
summary(m1, best1$posMinLambda)
summary(m1, best2$posMinLambda)
par(mfrow = c(1,3)); plot(m1, xvar="lambda");
plot(m1, xvar="objective"); plot(m1, xvar="grad")
set.seed(1234)
n <- 1000
q <- 6
k <- 4
# x <- runif(n)
X <- matrix(abs(rnorm(n*q)), n, q)
rownames(X) <- 1:n
colnames(X) <- paste0("X", 1:q)
theta <- matrix(c(1, 2, 0, 0, 0,
2, 0, 1, 0, 0,
0, 0, 0, 0, 0,
1, 0, 0, 1, -1.2,
0, 0, 0, 0, 0,
1.5, 0, .5, 0, 0,
0, 0, 0, 0, 0),
ncol=(k+1), byrow=TRUE)
rownames(theta) <- c("(intercept)", paste0("X", 1:q))
colnames(theta) <- c("(intercept)", "qnorm(p)", "p", "log(p)", "log(1-p)")
B <- function(p, k){matrix(cbind(1, qnorm(p), p, log(p), log(1-p)), nrow=(k+1), byrow=TRUE)}
Q <- function(p, theta, B, k, X){rowSums(X * t(theta %*% B(p, k)))}
pp <- runif(n)
y <- Q(p=pp, theta=theta, B=B, k=k, X=cbind(1, X))
s <- matrix(1, q+1, k+1); s[2:(q+1), 2] <- 0; s[1, 3:(k+1)] <- 0; s[2:3, 4:5] <- 0
s[4:5, 3] <- 0; s[6:7, 4:5] <- 0
m2 <- piqr(y ~ X, formula.p = ~ qnorm(p) + p + I(log(p)) + I(log(1-p)), s=s)
best1 <- gof.piqr(m2, method="AIC", plot=FALSE)
best2 <- gof.piqr(m2, method="BIC", plot=FALSE)
summary(m2, best1$posMinLambda)
summary(m2, best2$posMinLambda)
par(mfrow = c(1,3)); plot(m2, xvar="lambda");
plot(m2, xvar="objective"); plot(m2, xvar="grad")
## End(Not run)
# see the documentation for 'summary.piqr', and 'plot.piqr'
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.