Title: | Parametric Quantile Regression Models for Bounded Data |
---|---|
Description: | A collection of parametric quantile regression models for bounded data. At present, the package provides 13 parametric quantile regression models. It can specify regression structure for any quantile and shape parameters. It also provides several S3 methods to extract information from fitted model, such as residual analysis, prediction, plotting, and model comparison. For more computation efficient the [dpqr]'s, likelihood, score and hessian functions are written in C++. For further details see Mazucheli et. al (2022) <doi:10.1016/j.cmpb.2022.106816>. |
Authors: | André F. B. Menezes [aut, cre]
|
Maintainer: | André F. B. Menezes <[email protected]> |
License: | Apache License (>= 2) |
Version: | 0.0.6 |
Built: | 2025-03-05 04:49:42 UTC |
Source: | https://github.com/andrmenezes/unitquantreg |
The unitquantreg R package provides a collection of parametric quantile regression models for bounded data. At present, the package provides 13 parametric quantile regression models. It also enables several S3 methods to extract information from fitted model, such as residual analysis, prediction, plotting, and model comparison.
André F. B. Menezes [email protected]
Josmar Mazucheli [email protected]
Density function, distribution function, quantile function and random number generation function
for the arcsecant hyperbolic Weibull distribution reparametrized in terms of the -th quantile,
.
dashw(x, mu, theta, tau = 0.5, log = FALSE) pashw(q, mu, theta, tau = 0.5, lower.tail = TRUE, log.p = FALSE) qashw(p, mu, theta, tau = 0.5, lower.tail = TRUE, log.p = FALSE) rashw(n, mu, theta, tau = 0.5)
dashw(x, mu, theta, tau = 0.5, log = FALSE) pashw(q, mu, theta, tau = 0.5, lower.tail = TRUE, log.p = FALSE) qashw(p, mu, theta, tau = 0.5, lower.tail = TRUE, log.p = FALSE) rashw(n, mu, theta, tau = 0.5)
x , q
|
vector of positive quantiles. |
mu |
location parameter indicating the |
theta |
shape parameter. |
tau |
the parameter to specify which quantile use in the parametrization. |
log , log.p
|
logical; If TRUE, probabilities p are given as log(p). |
lower.tail |
logical; If TRUE, (default), |
p |
vector of probabilities. |
n |
number of observations. If |
Probability density function
Cumulative distribution function
Quantile function
Reparameterization
where is the shape parameter and
.
dashw
gives the density, pashw
gives the distribution function,
qashw
gives the quantile function and rashw
generates random deviates.
Invalid arguments will return an error message.
Josmar Mazucheli
André F. B. Menezes
Korkmaz, M. C., Chesneau, C. and Korkmaz, Z. S., (2021). A new alternative quantile regression model for the bounded response with educational measurements applications of OECD countries. Journal of Applied Statistics, 1–25.
set.seed(6969) x <- rashw(n = 1000, mu = 0.5, theta = 2.5, tau = 0.5) R <- range(x) S <- seq(from = R[1L], to = R[2L], by = 0.01) hist(x, prob = TRUE, main = 'arcsecant hyperbolic Weibull') lines(S, dashw(x = S, mu = 0.5, theta = 2.5, tau = 0.5), col = 2) plot(ecdf(x)) lines(S, pashw(q = S, mu = 0.5, theta = 2.5, tau = 0.5), col = 2) plot(quantile(x, probs = S), type = "l") lines(qashw(p = S, mu = 0.5, theta = 2.5, tau = 0.5), col = 2)
set.seed(6969) x <- rashw(n = 1000, mu = 0.5, theta = 2.5, tau = 0.5) R <- range(x) S <- seq(from = R[1L], to = R[2L], by = 0.01) hist(x, prob = TRUE, main = 'arcsecant hyperbolic Weibull') lines(S, dashw(x = S, mu = 0.5, theta = 2.5, tau = 0.5), col = 2) plot(ecdf(x)) lines(S, pashw(q = S, mu = 0.5, theta = 2.5, tau = 0.5), col = 2) plot(quantile(x, probs = S), type = "l") lines(qashw(p = S, mu = 0.5, theta = 2.5, tau = 0.5), col = 2)
The body fat percentage of individuals assisted in a public hospital in Curitiba, Paraná, Brasil.
data(bodyfat, package = "unitquantreg")
data(bodyfat, package = "unitquantreg")
A data.frame
with 298 observations and 9 columns:
arms
: Arms fat percentage.
legs
: Legs fat percentage.
body
: Body fat percentage.
android
: Android fat percentage.
gynecoid
: Ginecoid fat percentage.
bmi
: Body mass index - 24.71577.
age
: Age - 46.00.
sex
: Sex of individual. Female or male.
ipaq
: Factor variable indicating the sedentary, insufficiently
active or active.
André F. B. Menezes
Josmar Mazucheli
http://www.leg.ufpr.br/doku.php/publications:papercompanions:multquasibeta
Petterle, R. R., Bonat, W. H., Scarpin, C. T., Jonasson, T., and Borba, V. Z. C., (2020). Multivariate quasi-beta regression models for continuous bounded data. The International Journal of Biostatistics, 1–15, (preprint).
Mazucheli, J., Leiva, V., Alves, B., and Menezes A. F. B., (2021). A new quantile regression for modeling bounded data under a unit Birnbaum-Saunders distribution with applications in medicine and politics. Symmetry, 13(4) 1–21.
unitquantreg
objectsProduces a (half-)normal probability plot from a fitted model
object of class unitquantreg
.
hnp(object, ...) ## S3 method for class 'unitquantreg' hnp( object, nsim = 99, halfnormal = TRUE, plot = TRUE, output = TRUE, level = 0.95, resid.type = c("quantile", "cox-snell"), ... )
hnp(object, ...) ## S3 method for class 'unitquantreg' hnp( object, nsim = 99, halfnormal = TRUE, plot = TRUE, output = TRUE, level = 0.95, resid.type = c("quantile", "cox-snell"), ... )
object |
fitted model object of class |
... |
currently not used. |
nsim |
number of simulations used to compute envelope. Default is 99. |
halfnormal |
logical. If |
plot |
Should the (half-)normal plot be plotted? Default is |
output |
Should the output be returned? Default is |
level |
confidence level of the simulated envelope. Default is 0.95. |
resid.type |
type of residuals to be used. The default is |
Residuals plots with simulated envelope were proposed by Atkinson (1981) and can be construct as follows:
generate sample set of independent observations from the estimated
parameters of the fitted model;
fit the model using the generated sample, if halfnormal
is
TRUE
compute the absolute values of the residuals and arrange them in order;
repeat steps (1) and (2) nsim
number of times;
consider the sets of the
nsim
ordered statistics
of the residuals, then for each set compute the quantile level
/2,
the median and the quantile 1 - level
/2;
plot these values and the ordered residuals of the original sample set versus the expected order statistics of a (half)-normal distribution, which is approximated as
for half-normal plots, i.e., halfnormal=TRUE
or
for normal plots, i.e., halfnormal=FALSE
, where is the the
cumulative distribution function of standard Normal distribution for
quantile
residuals or the standard exponential distribution for the
cox-snell
residuals.
According to Atkinson (1981), if the model was correctly specified then no
more than level
100% of the observations are expected to appear
outside the envelope bands. Additionally, if a large proportion of the
observations lies outside the envelope, thus one has evidence against
the adequacy of the fitted model.
A list with the following components in ordered
(and absolute if halfnormal
is TRUE
) values:
obs |
the observed residuals. |
teo |
the theoretical residuals. |
lower |
lower envelope band. |
median |
median envelope band. |
upper |
upper envelope band. |
time_elapsed |
time elapsed to fit the |
André F. B. Menezes
Atkinson, A. C., (1981). Two graphical displays for outlying and influential observations in regression. Biometrika 68(1), 13–20.
Density function, distribution function, quantile function and random number generation function
for the Johnson SB distribution reparametrized in terms of the -th quantile,
.
djohnsonsb(x, mu, theta, tau = 0.5, log = FALSE) pjohnsonsb(q, mu, theta, tau = 0.5, lower.tail = TRUE, log.p = FALSE) qjohnsonsb(p, mu, theta, tau = 0.5, lower.tail = TRUE, log.p = FALSE) rjohnsonsb(n, mu, theta, tau = 0.5)
djohnsonsb(x, mu, theta, tau = 0.5, log = FALSE) pjohnsonsb(q, mu, theta, tau = 0.5, lower.tail = TRUE, log.p = FALSE) qjohnsonsb(p, mu, theta, tau = 0.5, lower.tail = TRUE, log.p = FALSE) rjohnsonsb(n, mu, theta, tau = 0.5)
x , q
|
vector of positive quantiles. |
mu |
location parameter indicating the |
theta |
nonnegative shape parameter. |
tau |
the parameter to specify which quantile is to used. |
log , log.p
|
logical; If TRUE, probabilities p are given as log(p). |
lower.tail |
logical; If TRUE, (default), |
p |
vector of probabilities. |
n |
number of observations. If |
Probability density function
Cumulative distribution function
Quantile function
Reparameterization
djohnsonsb
gives the density, pjohnsonsb
gives the distribution function,
qjohnsonsb
gives the quantile function and rjohnsonsb
generates random deviates.
Invalid arguments will return an error message.
Josmar Mazucheli
André F. B. Menezes
Lemonte, A. J. and Bazán, J. L., (2015). New class of Johnson SB distributions and its associated regression model for rates and proportions. Biometrical Journal, 58(4), 727–746.
Johnson, N. L., (1949). Systems of frequency curves generated by methods of translation. Biometrika, 36(1), 149–176.
set.seed(123) x <- rjohnsonsb(n = 1000, mu = 0.5, theta = 1.5, tau = 0.5) R <- range(x) S <- seq(from = R[1], to = R[2], by = 0.01) hist(x, prob = TRUE, main = 'Johnson SB') lines(S, djohnsonsb(x = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2) plot(ecdf(x)) lines(S, pjohnsonsb(q = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2) plot(quantile(x, probs = S), type = "l") lines(qjohnsonsb(p = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2)
set.seed(123) x <- rjohnsonsb(n = 1000, mu = 0.5, theta = 1.5, tau = 0.5) R <- range(x) S <- seq(from = R[1], to = R[2], by = 0.01) hist(x, prob = TRUE, main = 'Johnson SB') lines(S, djohnsonsb(x = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2) plot(ecdf(x)) lines(S, pjohnsonsb(q = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2) plot(quantile(x, probs = S), type = "l") lines(qjohnsonsb(p = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2)
Density function, distribution function, quantile function and random number generation for the Kumaraswamy distribution reparametrized in terms of the -th quantile,
.
dkum(x, mu, theta, tau = 0.5, log = FALSE) pkum(q, mu, theta, tau = 0.5, lower.tail = TRUE, log.p = FALSE) qkum(p, mu, theta, tau = 0.5, lower.tail = TRUE, log.p = FALSE) rkum(n, mu, theta, tau = 0.5)
dkum(x, mu, theta, tau = 0.5, log = FALSE) pkum(q, mu, theta, tau = 0.5, lower.tail = TRUE, log.p = FALSE) qkum(p, mu, theta, tau = 0.5, lower.tail = TRUE, log.p = FALSE) rkum(n, mu, theta, tau = 0.5)
x , q
|
vector of positive quantiles. |
mu |
location parameter indicating the |
theta |
nonnegative shape parameter. |
tau |
the parameter to specify which quantile is to used. |
log , log.p
|
logical; If TRUE, probabilities p are given as log(p). |
lower.tail |
logical; If TRUE, (default), |
p |
vector of probabilities. |
n |
number of observations. If |
Probability density function
Cumulative distribution function
Quantile function
Reparameterization
dkum
gives the density, pkum
gives the distribution function,
qkum
gives the quantile function and rkum
generates random deviates.
Invalid arguments will return an error message.
Josmar Mazucheli
André F. B. Menezes
Kumaraswamy, P., (1980). A generalized probability density function for double-bounded random processes. Journal of Hydrology, 46(1), 79–88.
Jones, M. C., (2009). Kumaraswamy's distribution: A beta-type distribution with some tractability advantages. Statistical Methodology, 6(1), 70-81.
set.seed(123) x <- rkum(n = 1000, mu = 0.5, theta = 1.5, tau = 0.5) R <- range(x) S <- seq(from = R[1], to = R[2], by = 0.01) hist(x, prob = TRUE, main = 'Kumaraswamy') lines(S, dkum(x = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2) plot(ecdf(x)) lines(S, pkum(q = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2) plot(quantile(x, probs = S), type = "l") lines(qkum(p = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2)
set.seed(123) x <- rkum(n = 1000, mu = 0.5, theta = 1.5, tau = 0.5) R <- range(x) S <- seq(from = R[1], to = R[2], by = 0.01) hist(x, prob = TRUE, main = 'Kumaraswamy') lines(S, dkum(x = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2) plot(ecdf(x)) lines(S, pkum(q = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2) plot(quantile(x, probs = S), type = "l") lines(qkum(p = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2)
Density function, distribution function, quantile function and random number generation function
for the Log-extended exponential-geometric distribution reparametrized in terms of the -th quantile,
.
dleeg(x, mu, theta, tau = 0.5, log = FALSE) pleeg(q, mu, theta, tau = 0.5, lower.tail = TRUE, log.p = FALSE) qleeg(p, mu, theta, tau = 0.5, lower.tail = TRUE, log.p = FALSE) rleeg(n, mu, theta, tau = 0.5)
dleeg(x, mu, theta, tau = 0.5, log = FALSE) pleeg(q, mu, theta, tau = 0.5, lower.tail = TRUE, log.p = FALSE) qleeg(p, mu, theta, tau = 0.5, lower.tail = TRUE, log.p = FALSE) rleeg(n, mu, theta, tau = 0.5)
x , q
|
vector of positive quantiles. |
mu |
location parameter indicating the |
theta |
nonnegative shape parameter. |
tau |
the parameter to specify which quantile is to be used. |
log , log.p
|
logical; If TRUE, probabilities p are given as log(p). |
lower.tail |
logical; If TRUE, (default), |
p |
vector of probabilities. |
n |
number of observations. If |
Probability density function
Cumulative distribution function
Quantile function
Reparameterization
dleeg
gives the density, pleeg
gives the distribution function,
qleeg
gives the quantile function and rleeg
generates random deviates.
Invalid arguments will return an error message.
Josmar Mazucheli [email protected]
André F. B. Menezes [email protected]
Jodrá, P. and Jiménez-Gamero, M. D., (2020). A quantile regression model for bounded responses based on the exponential-geometric distribution. Revstat - Statistical Journal, 18(4), 415–436.
set.seed(123) x <- rleeg(n = 1000, mu = 0.5, theta = 1.5, tau = 0.5) R <- range(x) S <- seq(from = R[1], to = R[2], by = 0.01) hist(x, prob = TRUE, main = 'Log-extended exponential-geometric') lines(S, dleeg(x = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2) plot(ecdf(x)) lines(S, pleeg(q = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2) plot(quantile(x, probs = S), type = "l") lines(qleeg(p = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2)
set.seed(123) x <- rleeg(n = 1000, mu = 0.5, theta = 1.5, tau = 0.5) R <- range(x) S <- seq(from = R[1], to = R[2], by = 0.01) hist(x, prob = TRUE, main = 'Log-extended exponential-geometric') lines(S, dleeg(x = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2) plot(ecdf(x)) lines(S, pleeg(q = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2) plot(quantile(x, probs = S), type = "l") lines(qleeg(p = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2)
unitquantreg
objects.Computes the likelihood-based statistics (Neg2LogLike, AIC, BIC and HQIC)
from unitquantreg
objects.
likelihood_stats(..., lt = NULL) ## S3 method for class 'likelihood_stats' print(x, ...)
likelihood_stats(..., lt = NULL) ## S3 method for class 'likelihood_stats' print(x, ...)
... |
|
lt |
a list with one or more |
x |
object of class |
Neg2LogLike: The log-likelihood is reported as
AIC: The Akaike's information criterion (AIC) is defined as
BIC: The Schwarz Bayesian information criterion (BIC) is defined as
HQIC: The Hannan and Quinn information criterion (HQIC) is defined as
where is the likelihood function.
A list with class "likelihood_stats"
containing the following
components:
call |
the matched call. |
stats |
ordered matrix according AIC value containg the likelihood based statistics. |
André F. B. Menezes
Josmar Mazucheli
Akaike, H. (1974). A new look at the statistical model identification. IEEE Transaction on Automatic Control, 19(6), 716–723.
Hannan, E. J. and Quinn, B. G. (1979). The determination of the order of an autoregression. Journal of the Royal Statistical Society, Series B, 41(2), 190–195.
Schwarz, G. (1978). Estimating the dimension of a model. Annals of Statistics, 6(2), 461–464.
data(sim_bounded, package = "unitquantreg") sim_bounded_curr <- sim_bounded[sim_bounded$family == "uweibull", ] models <- c("uweibull", "kum", "ulogistic") lt_fits <- lapply(models, function(fam) { unitquantreg(formula = y1 ~ x, tau = 0.5, data = sim_bounded_curr, family = fam) }) ans <- likelihood_stats(lt = lt_fits) ans
data(sim_bounded, package = "unitquantreg") sim_bounded_curr <- sim_bounded[sim_bounded$family == "uweibull", ] models <- c("uweibull", "kum", "ulogistic") lt_fits <- lapply(models, function(fam) { unitquantreg(formula = y1 ~ x, tau = 0.5, data = sim_bounded_curr, family = fam) }) ans <- likelihood_stats(lt = lt_fits) ans
Internal functions using in unitquantreg.fit
to compute the negative log-likelihood function, the score vector and the hessian
matrix using analytic expressions written in C++
.
loglike_unitquantreg(par, tau, family, linkobj, linkobj.theta, X, Z, y)
loglike_unitquantreg(par, tau, family, linkobj, linkobj.theta, X, Z, y)
par |
vector of regression model coefficients for |
tau |
quantile level, value between 0 and 1. |
family |
specify the distribution family name. |
linkobj , linkobj.theta
|
a function, usually obtained from
|
X |
design matrix related to the |
Z |
design matrix related to the |
y |
vector of response variable. |
unitquantreg
and unitquantregs
objectsMethods for extracting information from fitted regression models
objects of class unitquantreg
and unitquantregs
.
## S3 method for class 'unitquantreg' print(x, digits = max(4, getOption("digits") - 3), ...) ## S3 method for class 'unitquantreg' summary(object, correlation = FALSE, ...) ## S3 method for class 'unitquantreg' coef(object, type = c("full", "quantile", "shape"), ...) ## S3 method for class 'unitquantreg' vcov(object, ...) ## S3 method for class 'unitquantreg' logLik(object, ...) ## S3 method for class 'unitquantreg' confint(object, parm, level = 0.95, ...) ## S3 method for class 'unitquantreg' fitted(object, type = c("quantile", "shape", "full"), ...) ## S3 method for class 'unitquantreg' terms(x, type = c("quantile", "shape"), ...) ## S3 method for class 'unitquantreg' model.frame(formula, ...) ## S3 method for class 'unitquantreg' model.matrix(object, type = c("quantile", "shape"), ...) ## S3 method for class 'unitquantreg' update(object, formula., ..., evaluate = TRUE) ## S3 method for class 'unitquantregs' print(x, digits = max(3, getOption("digits") - 3), ...) ## S3 method for class 'unitquantregs' summary(object, digits = max(3, getOption("digits") - 3), ...)
## S3 method for class 'unitquantreg' print(x, digits = max(4, getOption("digits") - 3), ...) ## S3 method for class 'unitquantreg' summary(object, correlation = FALSE, ...) ## S3 method for class 'unitquantreg' coef(object, type = c("full", "quantile", "shape"), ...) ## S3 method for class 'unitquantreg' vcov(object, ...) ## S3 method for class 'unitquantreg' logLik(object, ...) ## S3 method for class 'unitquantreg' confint(object, parm, level = 0.95, ...) ## S3 method for class 'unitquantreg' fitted(object, type = c("quantile", "shape", "full"), ...) ## S3 method for class 'unitquantreg' terms(x, type = c("quantile", "shape"), ...) ## S3 method for class 'unitquantreg' model.frame(formula, ...) ## S3 method for class 'unitquantreg' model.matrix(object, type = c("quantile", "shape"), ...) ## S3 method for class 'unitquantreg' update(object, formula., ..., evaluate = TRUE) ## S3 method for class 'unitquantregs' print(x, digits = max(3, getOption("digits") - 3), ...) ## S3 method for class 'unitquantregs' summary(object, digits = max(3, getOption("digits") - 3), ...)
digits |
minimal number of significant digits. |
... |
additional argument(s) for methods. Currently not used. |
object , x
|
fitted model object of class |
correlation |
logical; if |
type |
character indicating type of fitted values to return. |
parm |
a specification of which parameters are to be given confidence intervals, either a vector of numbers or a vector of names. If missing, all parameters are considered. |
level |
the confidence level required. |
formula |
an R formula. |
formula. |
Changes to the formula see |
evaluate |
If true evaluate the new call else return the call. |
The summary
method gives Wald tests for the regressions coefficients
based on the observed Fisher information matrix. As usual the summary
method returns a list with relevant model statistics and estimates, which
can be printed using the print
method.
The coef
, vcov
, confint
and fitted
methods can
be use to extract, respectively, the estimated coefficients, the
estimated covariance matrix, the Wald confidence intervals, and fitted
values.
A logLik
method is also provide, then the AIC
function can be use to calculated the Akaike Information Criterion.
The generic methods terms
, model.frame
,
model.matrix
, update
and are also provided.
André F. B. Menezes
data(sim_bounded, package = "unitquantreg") sim_bounded_curr <- sim_bounded[sim_bounded$family == "uweibull", ] fit_1 <- unitquantreg(formula = y1 ~ x + z + I(x^2) | z + x, data = sim_bounded_curr, family = "uweibull", tau = 0.5, link.theta = "log") fit_1 summary(fit_1) vcov(fit_1) coef(fit_1) confint(fit_1) terms(fit_1) model.frame(fit_1)[1:5, ] model.matrix(fit_1)[1:5, ] update(fit_1, . ~ . -x) update(fit_1, . ~ . -z) update(fit_1, . ~ . -I(x^2)) update(fit_1, . ~ . | . -z) update(fit_1, . ~ . | . -x)
data(sim_bounded, package = "unitquantreg") sim_bounded_curr <- sim_bounded[sim_bounded$family == "uweibull", ] fit_1 <- unitquantreg(formula = y1 ~ x + z + I(x^2) | z + x, data = sim_bounded_curr, family = "uweibull", tau = 0.5, link.theta = "log") fit_1 summary(fit_1) vcov(fit_1) coef(fit_1) confint(fit_1) terms(fit_1) model.frame(fit_1)[1:5, ] model.matrix(fit_1)[1:5, ] update(fit_1, . ~ . -x) update(fit_1, . ~ . -z) update(fit_1, . ~ . -I(x^2)) update(fit_1, . ~ . | . -z) update(fit_1, . ~ . | . -x)
Calculate pairwise comparisons between fitted models performing
vuong test for objects of class unitquantreg
.
pairwise.vuong.test( ..., lt, p.adjust.method = p.adjust.methods, alternative = c("two.sided", "less", "greater") )
pairwise.vuong.test( ..., lt, p.adjust.method = p.adjust.methods, alternative = c("two.sided", "less", "greater") )
... |
|
lt |
a list with one or more |
p.adjust.method |
a character string specifying the method for multiple
testing adjustment; almost always one of
|
alternative |
indicates the alternative hypothesis and must be one
of |
Object of class "pairwise.htest"
data(sim_bounded, package = "unitquantreg") sim_bounded_curr <- sim_bounded[sim_bounded$family == "uweibull", ] models <- c("uweibull", "kum", "ulogistic") lt_fits <- lapply(models, function(fam) { unitquantreg(formula = y1 ~ x, tau = 0.5, data = sim_bounded_curr, family = fam) }) ans <- pairwise.vuong.test(lt = lt_fits) ans
data(sim_bounded, package = "unitquantreg") sim_bounded_curr <- sim_bounded[sim_bounded$family == "uweibull", ] models <- c("uweibull", "kum", "ulogistic") lt_fits <- lapply(models, function(fam) { unitquantreg(formula = y1 ~ x, tau = 0.5, data = sim_bounded_curr, family = fam) }) ans <- pairwise.vuong.test(lt = lt_fits) ans
unitquantreg
objectsProvide diagnostic plots to check model assumptions for fitted model
of class unitquantreg
.
## S3 method for class 'unitquantreg' plot( x, which = 1L:4L, caption = c("Residuals vs. indices of obs.", "Residuals vs. linear predictor", "Working response vs. linear predictor", "Half-normal plot of residuals"), sub.caption = paste(deparse(x$call), collapse = "\n"), main = "", ask = prod(par("mfcol")) < length(which) && dev.interactive(), ..., add.smooth = getOption("add.smooth"), type = "quantile", nsim = 99L, level = 0.95 )
## S3 method for class 'unitquantreg' plot( x, which = 1L:4L, caption = c("Residuals vs. indices of obs.", "Residuals vs. linear predictor", "Working response vs. linear predictor", "Half-normal plot of residuals"), sub.caption = paste(deparse(x$call), collapse = "\n"), main = "", ask = prod(par("mfcol")) < length(which) && dev.interactive(), ..., add.smooth = getOption("add.smooth"), type = "quantile", nsim = 99L, level = 0.95 )
x |
fitted model object of class |
which |
integer. if a subset of the plots is required, specify a subset of the numbers 1 to 4, see below for further details. |
caption |
character. Captions to appear above the plots. |
sub.caption |
character. Common title-above figures if there are multiple. |
main |
character. Title to each plot in addition to the above caption. |
ask |
logical. If |
... |
other parameters to be passed through to plotting functions. |
add.smooth |
logical. Indicates if a smoother should be added to most plots |
type |
character. Indicates type of residual to be used, see
|
nsim |
integer. Number of simulations in half-normal plots, see
|
level |
numeric. Confidence level of the simulated envelope, see
|
The plot
method for unitquantreg
objects produces four types
of diagnostic plot.
The which
argument can be used to select a subset of currently four
supported plot, which are: Residuals versus indices of observations
(which = 1
); Residuals versus linear predictor (which = 2
);
Working response versus linear predictor (which = 3
) to
check possible misspecification of link function; Half-normal plot of
residuals (which = 4
) to check distribution assumption.
No return value, called for side effects.
André F. B. Menezes
Dunn, P. K. and Smyth, G. K. (2018) Generalized Linear Models With Examples in R, Springer, New York.
residuals.unitquantreg
,
hnp.unitquantreg
,
unitquantreg
.
unitquantregs
objectsProvide two type of plots for unitquantregs
objects.
## S3 method for class 'unitquantregs' plot( x, which = c("coef", "conddist"), output_df = FALSE, parm = NULL, level = 0.95, mean_effect = FALSE, mfrow = NULL, mar = NULL, ylim = NULL, main = NULL, col = gray(c(0, 0.75)), border = NULL, cex = 1, pch = 20, type = "b", xlab = bquote("Quantile level (" * tau * ")"), ylab = "Estimate effect", dist_type = c("density", "cdf"), at_avg = TRUE, at_obs = NULL, legend_position = "topleft", ... )
## S3 method for class 'unitquantregs' plot( x, which = c("coef", "conddist"), output_df = FALSE, parm = NULL, level = 0.95, mean_effect = FALSE, mfrow = NULL, mar = NULL, ylim = NULL, main = NULL, col = gray(c(0, 0.75)), border = NULL, cex = 1, pch = 20, type = "b", xlab = bquote("Quantile level (" * tau * ")"), ylab = "Estimate effect", dist_type = c("density", "cdf"), at_avg = TRUE, at_obs = NULL, legend_position = "topleft", ... )
x |
fitted model object of class |
which |
character. Indicate the type of plot. Currently supported are |
output_df |
logical. Should |
parm |
a specification of which parameters are to be plotted, either a vector of numbers or a vector of names. By default, all parameters are considered. |
level |
level of significance for the confidence interval of parameters. |
mean_effect |
logical. Should a line for the mean effect coefficients be added? |
mfrow , mar , ylim , main , col , border , cex , pch , type , xlab , ylab
|
graphical parameters. |
dist_type |
character. Which conditional distribution should be plotted?
The options are |
at_avg |
logical. Should consider the conditional distribution at average values of covariates? |
at_obs |
list. List with name and values for each covariate. |
legend_position |
character. The legend position argument used in |
... |
other parameters to be passed through to plotting functions. |
The plot method for unitquantregs
objects is inspired in PROC QUANTREG of SAS/STAT.
This plot method provide two type of visualizations.
If which = "coef"
plot the estimated coefficients for several quantiles.
If which = "conddist"
plot the conditional distribution at specific values of
covariates. The conditional distribution could be the cumulative distribution function
if dist_type = "cdf"
or the probability density function if dist_type = "pdf"
.
If output_df = TRUE
then returns a data.frame used to plot.
Otherwise, no return value, called for side effects.
André F. B. Menezes
unitquantreg
classExtract various types of predictions from unit quantile regression models.
## S3 method for class 'unitquantreg' predict( object, newdata, type = c("link", "quantile", "shape", "terms"), interval = c("none", "confidence"), level = 0.95, se.fit = FALSE, ... )
## S3 method for class 'unitquantreg' predict( object, newdata, type = c("link", "quantile", "shape", "terms"), interval = c("none", "confidence"), level = 0.95, se.fit = FALSE, ... )
object |
fitted model object of class |
newdata |
optionally, a data frame in which to look for variables with which to predict. If omitted, the original observations are used. |
type |
character indicating type of predictions. The options are
|
interval |
type of interval desired. The options are |
level |
coverage probability for the confidence intervals. Default is
|
se.fit |
logical. If |
... |
currently not used. |
If se.fit = FALSE
then returns a data.frame
with
predict values and confidence interval if interval = TRUE
.
If se.fit = TRUE
returns a list with components:
fit |
Predictions, as for |
se.fit |
Estimated standard errors. |
For type = "terms"
the output is a data.frame
with a columns
per term.
André F. B. Menezes
unitquantreg
objectsExtract various types of residuals from unit quantile regression models.
## S3 method for class 'unitquantreg' residuals(object, type = c("quantile", "cox-snell", "working", "partial"), ...)
## S3 method for class 'unitquantreg' residuals(object, type = c("quantile", "cox-snell", "working", "partial"), ...)
object |
fitted model object of class |
type |
character indicating type of residuals. The options are
|
... |
currently not used. |
The residuals
method can compute quantile
and Cox-Snell residuals. These residuals are defined, respectively, by
and
where and
are the fitted values
of parameters
and
,
is
the cumulative distribution function (c.d.f.) and
is the
c.d.f. of standard Normal distribution.
Apart from the variability due the estimates of parameters,if the fitted
regression model is correctly specified then the quantile
residuals, , follow a standard Normal distribution and
the Cox-Snell residuals,
, follow a standard exponential
distribution.
Numeric vector of residuals extract from an object of class
unitquantreg
.
André F. B. Menezes
Cox, D. R. and Snell E. J., (1968). A general definition of residuals. Journal of the Royal Statistical Society - Series B, 30(2), 248–265.
Dunn, P. K. and Smyth, G. K., (1996). Randomized quantile residuals. Journal of Computational and Graphical Statistics, 5(3), 236–244.
This data set was simulated from all families of distributions
available in unitquantreg
package considering
the median, i.e., .
data(sim_bounded, package = "unitquantreg")
data(sim_bounded, package = "unitquantreg")
data.frame
with 1300 observations and 5 columns:
y1
: simulated response variable with constant shape parameter,
.
y2
: simulated response variable with regression structure in
the shape parameter, ), where
.
x
: covariate related to , i.e., the median.
z
: covariate related to , i.e., the shape parameter.
family
: string indicating the family of distribution.
There are two response variable, namely y1
and y2
.
The former was simulated considering a regression structure for
and one covariate simulated from a standard uniform distribution,
where the true vector of coefficients for
is
and
.
The latter was simulated assuming a regression structure for both
and
(shape parameter) and only one independent covariates
simulated from two standard uniform distributions. The true vectors of
coefficients for
and
are
and
,
respectively.
André F. B. Menezes
Density function, distribution function, quantile function and random number generation function
for the unit-Birnbaum-Saunders distribution reparametrized in terms of the -th quantile,
.
dubs(x, mu, theta, tau = 0.5, log = FALSE) pubs(q, mu, theta, tau = 0.5, lower.tail = TRUE, log.p = FALSE) qubs(p, mu, theta, tau = 0.5, lower.tail = TRUE, log.p = FALSE) rubs(n, mu, theta, tau = 0.5)
dubs(x, mu, theta, tau = 0.5, log = FALSE) pubs(q, mu, theta, tau = 0.5, lower.tail = TRUE, log.p = FALSE) qubs(p, mu, theta, tau = 0.5, lower.tail = TRUE, log.p = FALSE) rubs(n, mu, theta, tau = 0.5)
x , q
|
vector of positive quantiles. |
mu |
location parameter indicating the |
theta |
nonnegative shape parameter. |
tau |
the parameter to specify which quantile is to be used. |
log , log.p
|
logical; If TRUE, probabilities p are given as log(p). |
lower.tail |
logical; If TRUE, (default), |
p |
vector of probabilities. |
n |
number of observations. If |
Probability density function
Cumulative distribution function
Quantile function
Reparameterization
where
dubs
gives the density, pubs
gives the distribution function,
qubs
gives the quantile function and rubs
generates random deviates.
Invalid arguments will return an error message.
Josmar Mazucheli [email protected]
André F. B. Menezes [email protected]
Birnbaum, Z. W. and Saunders, S. C., (1969). A new family of life distributions. Journal of Applied Probability, 6(2), 637–652. Mazucheli, J., Menezes, A. F. B. and Dey, S., (2018). The unit-Birnbaum-Saunders distribution with applications. Chilean Journal of Statistics, 9(1), 47–57.
Mazucheli, J., Alves, B. and Menezes, A. F. B., (2021). A new quantile regression for modeling bounded data under a unit Birnbaum-Saunders distribution with applications. Simmetry, (), 1–28.
set.seed(123) x <- rubs(n = 1000, mu = 0.5, theta = 1.5, tau = 0.5) R <- range(x) S <- seq(from = R[1], to = R[2], by = 0.01) hist(x, prob = TRUE, main = 'unit-Birnbaum-Saunders') lines(S, dubs(x = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2) plot(ecdf(x)) lines(S, pubs(q = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2) plot(quantile(x, probs = S), type = "l") lines(qubs(p = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2)
set.seed(123) x <- rubs(n = 1000, mu = 0.5, theta = 1.5, tau = 0.5) R <- range(x) S <- seq(from = R[1], to = R[2], by = 0.01) hist(x, prob = TRUE, main = 'unit-Birnbaum-Saunders') lines(S, dubs(x = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2) plot(ecdf(x)) lines(S, pubs(q = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2) plot(quantile(x, probs = S), type = "l") lines(qubs(p = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2)
Density function, distribution function, quantile function and random number generation function
for the unit-Burr-XII distribution reparametrized in terms of the -th quantile,
.
duburrxii(x, mu, theta, tau = 0.5, log = FALSE) puburrxii(q, mu, theta, tau = 0.5, lower.tail = TRUE, log.p = FALSE) quburrxii(p, mu, theta, tau = 0.5, lower.tail = TRUE, log.p = FALSE) ruburrxii(n, mu, theta, tau = 0.5)
duburrxii(x, mu, theta, tau = 0.5, log = FALSE) puburrxii(q, mu, theta, tau = 0.5, lower.tail = TRUE, log.p = FALSE) quburrxii(p, mu, theta, tau = 0.5, lower.tail = TRUE, log.p = FALSE) ruburrxii(n, mu, theta, tau = 0.5)
x , q
|
vector of positive quantiles. |
mu |
location parameter indicating the |
theta |
nonnegative shape parameter. |
tau |
the parameter to specify which quantile is to used. |
log , log.p
|
logical; If TRUE, probabilities p are given as log(p). |
lower.tail |
logical; If TRUE, (default), |
p |
vector of probabilities. |
n |
number of observations. If |
Probability density function
Cumulative distribution function
Quantile function
Reparameterization
duburrxii
gives the density, puburrxii
gives the distribution function,
quburrxii
gives the quantile function and ruburrxii
generates random deviates.
Invalid arguments will return an error message.
Josmar Mazucheli [email protected]
André F. B. Menezes [email protected]
Korkmaz M. C. and Chesneau, C., (2021). On the unit Burr-XII distribution with the quantile regression modeling and applications. Computational and Applied Mathematics, 40(29), 1–26.
set.seed(123) x <- ruburrxii(n = 1000, mu = 0.5, theta = 1.5, tau = 0.5) R <- range(x) S <- seq(from = R[1], to = R[2], by = 0.01) hist(x, prob = TRUE, main = 'unit-Burr-XII') lines(S, duburrxii(x = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2) plot(ecdf(x)) lines(S, puburrxii(q = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2) plot(quantile(x, probs = S), type = "l") lines(quburrxii(p = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2)
set.seed(123) x <- ruburrxii(n = 1000, mu = 0.5, theta = 1.5, tau = 0.5) R <- range(x) S <- seq(from = R[1], to = R[2], by = 0.01) hist(x, prob = TRUE, main = 'unit-Burr-XII') lines(S, duburrxii(x = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2) plot(ecdf(x)) lines(S, puburrxii(q = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2) plot(quantile(x, probs = S), type = "l") lines(quburrxii(p = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2)
Density function, distribution function, quantile function and random number generation function
for the unit-Chen distribution reparametrized in terms of the -th quantile,
.
duchen(x, mu, theta, tau = 0.5, log = FALSE) puchen(q, mu, theta, tau = 0.5, lower.tail = TRUE, log.p = FALSE) quchen(p, mu, theta, tau = 0.5, lower.tail = TRUE, log.p = FALSE) ruchen(n, mu, theta, tau = 0.5)
duchen(x, mu, theta, tau = 0.5, log = FALSE) puchen(q, mu, theta, tau = 0.5, lower.tail = TRUE, log.p = FALSE) quchen(p, mu, theta, tau = 0.5, lower.tail = TRUE, log.p = FALSE) ruchen(n, mu, theta, tau = 0.5)
x , q
|
vector of positive quantiles. |
mu |
location parameter indicating the |
theta |
nonnegative shape parameter. |
tau |
the parameter to specify which quantile is to be used. |
log , log.p
|
logical; If TRUE, probabilities p are given as log(p). |
lower.tail |
logical; If TRUE, (default), |
p |
vector of probabilities. |
n |
number of observations. If |
Probability density function
Cumulative distribution function
Quantile function
Reparameterization
duchen
gives the density, puchen
gives the distribution function,
quchen
gives the quantile function and ruchen
generates random deviates.
Invalid arguments will return an error message.
Josmar Mazucheli [email protected]
André F. B. Menezes [email protected]
Korkmaz, M. C., Emrah, A., Chesneau, C. and Yousof, H. M., (2020). On the unit-Chen distribution with associated quantile regression and applications. Journal of Applied Statistics, 44(1) 1–22.
set.seed(123) x <- ruchen(n = 1000, mu = 0.5, theta = 1.5, tau = 0.5) R <- range(x) S <- seq(from = R[1], to = R[2], by = 0.01) hist(x, prob = TRUE, main = 'unit-Chen') lines(S, duchen(x = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2) plot(ecdf(x)) lines(S, puchen(q = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2) plot(quantile(x, probs = S), type = "l") lines(quchen(p = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2)
set.seed(123) x <- ruchen(n = 1000, mu = 0.5, theta = 1.5, tau = 0.5) R <- range(x) S <- seq(from = R[1], to = R[2], by = 0.01) hist(x, prob = TRUE, main = 'unit-Chen') lines(S, duchen(x = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2) plot(ecdf(x)) lines(S, puchen(q = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2) plot(quantile(x, probs = S), type = "l") lines(quchen(p = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2)
Density function, distribution function, quantile function and random number generation function
for the unit-Half-Normal-E distribution reparametrized in terms of the -th quantile,
.
dughne(x, mu, theta, tau = 0.5, log = FALSE) pughne(q, mu, theta, tau = 0.5, lower.tail = TRUE, log.p = FALSE) qughne(p, mu, theta, tau = 0.5, lower.tail = TRUE, log.p = FALSE) rughne(n, mu, theta, tau = 0.5)
dughne(x, mu, theta, tau = 0.5, log = FALSE) pughne(q, mu, theta, tau = 0.5, lower.tail = TRUE, log.p = FALSE) qughne(p, mu, theta, tau = 0.5, lower.tail = TRUE, log.p = FALSE) rughne(n, mu, theta, tau = 0.5)
x , q
|
vector of positive quantiles. |
mu |
location parameter indicating the |
theta |
nonnegative shape parameter. |
tau |
the parameter to specify which quantile is to be used. |
log , log.p
|
logical; If TRUE, probabilities p are given as log(p). |
lower.tail |
logical; If TRUE, (default), |
p |
vector of probabilities. |
n |
number of observations. If |
Probability density function
Cumulative distribution function
Quantile function
Reparameterization
dughne
gives the density, pughne
gives the distribution function,
qughne
gives the quantile function and rughne
generates random deviates.
Invalid arguments will return an error message.
Josmar Mazucheli [email protected]
André F. B. Menezes [email protected]
Korkmaz, M. C., (2020). The unit generalized half normal distribution: A new bounded distribution with inference and application. University Politehnica of Bucharest Scientific, 82(2), 133–140.
set.seed(123) x <- rughne(n = 1000, mu = 0.5, theta = 2, tau = 0.5) R <- range(x) S <- seq(from = R[1], to = R[2], by = 0.01) hist(x, prob = TRUE, main = 'unit-Half-Normal-E') lines(S, dughne(x = S, mu = 0.5, theta = 2, tau = 0.5), col = 2) plot(ecdf(x)) lines(S, pughne(q = S, mu = 0.5, theta = 2, tau = 0.5), col = 2) plot(quantile(x, probs = S), type = "l") lines(qughne(p = S, mu = 0.5, theta = 2, tau = 0.5), col = 2)
set.seed(123) x <- rughne(n = 1000, mu = 0.5, theta = 2, tau = 0.5) R <- range(x) S <- seq(from = R[1], to = R[2], by = 0.01) hist(x, prob = TRUE, main = 'unit-Half-Normal-E') lines(S, dughne(x = S, mu = 0.5, theta = 2, tau = 0.5), col = 2) plot(ecdf(x)) lines(S, pughne(q = S, mu = 0.5, theta = 2, tau = 0.5), col = 2) plot(quantile(x, probs = S), type = "l") lines(qughne(p = S, mu = 0.5, theta = 2, tau = 0.5), col = 2)
Density function, distribution function, quantile function and random number generation function
for the unit-Half-Normal-X distribution reparametrized in terms of the -th quantile,
.
dughnx(x, mu, theta, tau = 0.5, log = FALSE) pughnx(q, mu, theta, tau = 0.5, lower.tail = TRUE, log.p = FALSE) qughnx(p, mu, theta, tau = 0.5, lower.tail = TRUE, log.p = FALSE) rughnx(n, mu, theta, tau = 0.5)
dughnx(x, mu, theta, tau = 0.5, log = FALSE) pughnx(q, mu, theta, tau = 0.5, lower.tail = TRUE, log.p = FALSE) qughnx(p, mu, theta, tau = 0.5, lower.tail = TRUE, log.p = FALSE) rughnx(n, mu, theta, tau = 0.5)
x , q
|
vector of positive quantiles. |
mu |
location parameter indicating the |
theta |
nonnegative shape parameter. |
tau |
the parameter to specify which quantile is to be used. |
log , log.p
|
logical; If TRUE, probabilities p are given as log(p). |
lower.tail |
logical; If TRUE, (default), |
p |
vector of probabilities. |
n |
number of observations. If |
Probability density function
Cumulative density function
Quantile Function
Reparametrization
dughnx
gives the density, pughnx
gives the distribution function,
qughnx
gives the quantile function and rughnx
generates random deviates.
Invalid arguments will return an error message.
Josmar Mazucheli [email protected]
André F. B. Menezes [email protected]
Bakouch, H. S., Nik, A. S., Asgharzadeh, A. and Salinas, H. S., (2021). A flexible probability model for proportion data: Unit-Half-Normal distribution. Communications in Statistics: CaseStudies, Data Analysis and Applications, 0(0), 1–18.
set.seed(123) x <- rughnx(n = 1000, mu = 0.5, theta = 2, tau = 0.5) R <- range(x) S <- seq(from = R[1], to = R[2], by = 0.01) hist(x, prob = TRUE, main = 'unit-Half-Normal-X') lines(S, dughnx(x = S, mu = 0.5, theta = 2, tau = 0.5), col = 2) plot(ecdf(x)) lines(S, pughnx(q = S, mu = 0.5, theta = 2, tau = 0.5), col = 2) plot(quantile(x, probs = S), type = "l") lines(qughnx(p = S, mu = 0.5, theta = 2, tau = 0.5), col = 2)
set.seed(123) x <- rughnx(n = 1000, mu = 0.5, theta = 2, tau = 0.5) R <- range(x) S <- seq(from = R[1], to = R[2], by = 0.01) hist(x, prob = TRUE, main = 'unit-Half-Normal-X') lines(S, dughnx(x = S, mu = 0.5, theta = 2, tau = 0.5), col = 2) plot(ecdf(x)) lines(S, pughnx(q = S, mu = 0.5, theta = 2, tau = 0.5), col = 2) plot(quantile(x, probs = S), type = "l") lines(qughnx(p = S, mu = 0.5, theta = 2, tau = 0.5), col = 2)
Density function, distribution function, quantile function and random number deviates
for the unit-Gompertz distribution reparametrized in terms of the -th quantile,
.
dugompertz(x, mu, theta, tau = 0.5, log = FALSE) pugompertz(q, mu, theta, tau = 0.5, lower.tail = TRUE, log.p = FALSE) qugompertz(p, mu, theta, tau = 0.5, lower.tail = TRUE, log.p = FALSE) rugompertz(n, mu, theta, tau = 0.5)
dugompertz(x, mu, theta, tau = 0.5, log = FALSE) pugompertz(q, mu, theta, tau = 0.5, lower.tail = TRUE, log.p = FALSE) qugompertz(p, mu, theta, tau = 0.5, lower.tail = TRUE, log.p = FALSE) rugompertz(n, mu, theta, tau = 0.5)
x , q
|
vector of positive quantiles. |
mu |
location parameter indicating the |
theta |
nonnegative shape parameter. |
tau |
the parameter to specify which quantile is to be used. |
log , log.p
|
logical; If TRUE, probabilities p are given as log(p). |
lower.tail |
logical; If TRUE, (default), |
p |
vector of probabilities. |
n |
number of observations. If |
Probability density function
Cumulative density function
Quantile Function
Reparameterization
dugompertz
gives the density, pugompertz
gives the distribution function,
qugompertz
gives the quantile function and rugompertz
generates random deviates.
Invalid arguments will return an error message.
Josmar Mazucheli [email protected]
André F. B. Menezes [email protected]
Mazucheli, J., Menezes, A. F. and Dey, S., (2019). Unit-Gompertz Distribution with Applications. Statistica, 79(1), 25-43.
set.seed(123) x <- rugompertz(n = 1000, mu = 0.5, theta = 2, tau = 0.5) R <- range(x) S <- seq(from = R[1], to = R[2], by = 0.01) hist(x, prob = TRUE, main = 'unit-Gompertz') lines(S, dugompertz(x = S, mu = 0.5, theta = 2, tau = 0.5), col = 2) plot(ecdf(x)) lines(S, pugompertz(q = S, mu = 0.5, theta = 2, tau = 0.5), col = 2) plot(quantile(x, probs = S), type = "l") lines(qugompertz(p = S, mu = 0.5, theta = 2, tau = 0.5), col = 2)
set.seed(123) x <- rugompertz(n = 1000, mu = 0.5, theta = 2, tau = 0.5) R <- range(x) S <- seq(from = R[1], to = R[2], by = 0.01) hist(x, prob = TRUE, main = 'unit-Gompertz') lines(S, dugompertz(x = S, mu = 0.5, theta = 2, tau = 0.5), col = 2) plot(ecdf(x)) lines(S, pugompertz(q = S, mu = 0.5, theta = 2, tau = 0.5), col = 2) plot(quantile(x, probs = S), type = "l") lines(qugompertz(p = S, mu = 0.5, theta = 2, tau = 0.5), col = 2)
Density function, distribution function, quantile function and random number generation function
for the unit-Gumbel distribution reparametrized in terms of the -th quantile,
.
dugumbel(x, mu, theta, tau = 0.5, log = FALSE) pugumbel(q, mu, theta, tau = 0.5, lower.tail = TRUE, log.p = FALSE) qugumbel(p, mu, theta, tau = 0.5, lower.tail = TRUE, log.p = FALSE) rugumbel(n, mu, theta, tau = 0.5)
dugumbel(x, mu, theta, tau = 0.5, log = FALSE) pugumbel(q, mu, theta, tau = 0.5, lower.tail = TRUE, log.p = FALSE) qugumbel(p, mu, theta, tau = 0.5, lower.tail = TRUE, log.p = FALSE) rugumbel(n, mu, theta, tau = 0.5)
x , q
|
vector of positive quantiles. |
mu |
location parameter indicating the |
theta |
nonnegative shape parameter. |
tau |
the parameter to specify which quantile use in the parametrization. |
log , log.p
|
logical; If TRUE, probabilities p are given as log(p). |
lower.tail |
logical; If TRUE, (default), |
p |
vector of probabilities. |
n |
number of observations. If |
Probability density function
Cumulative distribution function
Quantile function
Reparameterization
where and
is the shape parameter.
dugumbel
gives the density, pugumbel
gives the distribution function,
qugumbel
gives the quantile function and rugumbel
generates random deviates.
Invalid arguments will return an error message.
Josmar Mazucheli
Andre F. B. Menezes
Mazucheli, J. and Alves, B., (2021). The unit-Gumbel Quantile Regression Model for Proportion Data. Under Review.
Gumbel, E. J., (1941). The return period of flood flows. The Annals of Mathematical Statistics, 12(2), 163–190.
set.seed(6969) x <- rugumbel(n = 1000, mu = 0.5, theta = 1.5, tau = 0.5) R <- range(x) S <- seq(from = R[1], to = R[2], by = 0.01) hist(x, prob = TRUE, main = 'unit-Gumbel') lines(S, dugumbel(x = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2) plot(ecdf(x)) lines(S, pugumbel(q = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2) plot(quantile(x, probs = S), type = "l") lines(qugumbel(p = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2)
set.seed(6969) x <- rugumbel(n = 1000, mu = 0.5, theta = 1.5, tau = 0.5) R <- range(x) S <- seq(from = R[1], to = R[2], by = 0.01) hist(x, prob = TRUE, main = 'unit-Gumbel') lines(S, dugumbel(x = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2) plot(ecdf(x)) lines(S, pugumbel(q = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2) plot(quantile(x, probs = S), type = "l") lines(qugumbel(p = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2)
Density function, distribution function, quantile function and random number generation for the unit-Logistic distribution reparametrized in terms of the -th quantile,
.
dulogistic(x, mu, theta, tau = 0.5, log = FALSE) pulogistic(q, mu, theta, tau = 0.5, lower.tail = TRUE, log.p = FALSE) qulogistic(p, mu, theta, tau = 0.5, lower.tail = TRUE, log.p = FALSE) rulogistic(n, mu, theta, tau = 0.5)
dulogistic(x, mu, theta, tau = 0.5, log = FALSE) pulogistic(q, mu, theta, tau = 0.5, lower.tail = TRUE, log.p = FALSE) qulogistic(p, mu, theta, tau = 0.5, lower.tail = TRUE, log.p = FALSE) rulogistic(n, mu, theta, tau = 0.5)
x , q
|
vector of positive quantiles. |
mu |
location parameter indicating the |
theta |
nonnegative shape parameter. |
tau |
the parameter to specify which quantile is to used. |
log , log.p
|
logical; If TRUE, probabilities p are given as log(p). |
lower.tail |
logical; If TRUE, (default), |
p |
vector of probabilities. |
n |
number of observations. If |
Probability density function
Cumulative distribution function
Quantile function
Reparameterization
dulogistic
gives the density, pulogistic
gives the distribution function,
qulogistic
gives the quantile function and rulogistic
generates random deviates.
Invalid arguments will return an error message.
Josmar Mazucheli [email protected]
André F. B. Menezes [email protected]
Paz, R. F., Balakrishnan, N. and Bazán, J. L., 2019. L-Logistic regression models: Prior sensitivity analysis, robustness to outliers and applications. Brazilian Journal of Probability and Statistics, 33(3), 455–479.
set.seed(123) x <- rulogistic(n = 1000, mu = 0.5, theta = 1.5, tau = 0.5) R <- range(x) S <- seq(from = R[1], to = R[2], by = 0.01) hist(x, prob = TRUE, main = 'unit-Logistic') lines(S, dulogistic(x = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2) plot(ecdf(x)) lines(S, pulogistic(q = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2) plot(quantile(x, probs = S), type = "l") lines(qulogistic(p = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2)
set.seed(123) x <- rulogistic(n = 1000, mu = 0.5, theta = 1.5, tau = 0.5) R <- range(x) S <- seq(from = R[1], to = R[2], by = 0.01) hist(x, prob = TRUE, main = 'unit-Logistic') lines(S, dulogistic(x = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2) plot(ecdf(x)) lines(S, pulogistic(q = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2) plot(quantile(x, probs = S), type = "l") lines(qulogistic(p = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2)
Fit a collection of parametric unit quantile regression model
by maximum likelihood using the log-likelihood function, the score vector
and the hessian matrix implemented in C++
.
unitquantreg( formula, data, subset, na.action, tau, family, link = c("logit", "probit", "cloglog", "cauchit"), link.theta = c("identity", "log", "sqrt"), start = NULL, control = unitquantreg.control(), model = TRUE, x = FALSE, y = TRUE ) unitquantreg.fit( y, X, Z = NULL, tau, family, link, link.theta, start = NULL, control = unitquantreg.control() )
unitquantreg( formula, data, subset, na.action, tau, family, link = c("logit", "probit", "cloglog", "cauchit"), link.theta = c("identity", "log", "sqrt"), start = NULL, control = unitquantreg.control(), model = TRUE, x = FALSE, y = TRUE ) unitquantreg.fit( y, X, Z = NULL, tau, family, link, link.theta, start = NULL, control = unitquantreg.control() )
formula |
symbolic description of the quantile model like |
data |
data.frame contain the variables in the model. |
subset |
an optional vector specifying a subset of observations to be used in the fitting process. |
na.action |
a function which indicates what should happen when the
data contain |
tau |
numeric vector. The quantile(s) to be estimated, i.e.,
number between 0 and 1. If just one quantile is specified an object of class
|
family |
character. Specify the distribution family. |
link |
character. Specify the link function in the quantile model.
Currently supported are |
link.theta |
character. Specify the link function in the shape model.
Currently supported are |
start |
numeric vector. An optional vector with starting values for all parameters. |
control |
list. Control arguments specified via |
model |
logical. Indicates whether model frame should be included as a component of the returned value. |
x , y
|
logical. If |
X , Z
|
numeric matrix. Regressor matrix for the quantile and shape model,
respectively. Default is constant shape model, i.e., |
The parameter estimation and inference are performed under the frequentist paradigm.
The optimx
R package is use, since allows different optimization
technique to maximize the log-likelihood function. The analytical score function are
use in the maximization and the standard errors are computed using the
analytical hessian matrix, both are implemented in efficient away using C++
.
unitquantreg
can return an object of
class unitquantreg
if tau
is a scalar, i.e., a list with
the following components.
family |
the distribution family name. |
coefficients |
a list with elements |
fitted.values |
a list with elements |
linear.predictors |
a list with elements |
link |
a list with elements |
tau |
the quantile specify. |
loglik |
log-likelihood of the fitted model. |
gradient |
gradient evaluate at maximum likelihood estimates. |
vcov |
covariance matrix of all parameters in the model. |
nobs |
number of observations. |
npar |
number of parameters. |
df.residual |
residual degrees of freedom in the fitted model. |
theta_const |
logical indicating if the |
control |
the control parameters used to fit the model. |
iterations |
number of iterations of optimization method. |
converged |
logical, if |
kkt |
a list of logical |
elapsed_time |
time elapsed to fit the model. |
call |
the original function call. |
formula |
the original model formula. |
terms |
a list with elements |
model |
the full model frame, if |
y |
the response vector, if |
x |
a list with elements |
While unitquantreg.fit
returns an unclassed list with
components up to elapsed_time
.
If tau
is a numeric vector with length greater than one an object of
class unitquantregs
is returned, which consist of list of objects of
class unitquantreg
for each specified quantiles.
André F. B. Menezes
Auxiliary function that control fitting of unit quantile
regression models using unitquantreg
.
unitquantreg.control( method = "BFGS", hessian = FALSE, gradient = TRUE, maxit = 5000, factr = 1e+07, reltol = sqrt(.Machine$double.eps), trace = 0L, starttests = FALSE, dowarn = FALSE, ... )
unitquantreg.control( method = "BFGS", hessian = FALSE, gradient = TRUE, maxit = 5000, factr = 1e+07, reltol = sqrt(.Machine$double.eps), trace = 0L, starttests = FALSE, dowarn = FALSE, ... )
method |
string. Specify the method argument passed to |
hessian |
logical. Should use the numerically Hessian matrix to compute
variance-covariance? Default is |
gradient |
logical. Should use the analytic gradient? Default is
|
maxit |
integer. Specify the maximal number of iterations passed to
|
factr |
numeric.Controls the convergence of the |
reltol |
numeric. Relative convergence tolerance passed to |
trace |
non-negative integer. If positive, tracing information on the progress of the optimization is produced. |
starttests |
logical. Should |
dowarn |
logical. Show warnings generated by |
... |
arguments passed to |
The control
argument of
unitquantreg
uses the arguments of
unitquantreg.control
. In particular, the
parameters in unitquantreg
are estimated by
maximum likelihood using the optimx
, which is a
general-purpose optimization wrapper function that calls other R tools for
optimization, including the existing optim
function.
The main advantage of optimx
is to unify the tools
allowing a number of different optimization methods and provide sanity checks.
A list with components named as the arguments.
André F. B. Menezes
Nash, J. C. and Varadhan, R. (2011). Unifying Optimization Algorithms to Aid Software System Users: optimx for R., Journal of Statistical Software, 43(9), 1–14.
optimx
for more details about control
parameters and unitquantreg.fit
the fitting
procedure used by unitquantreg
.
data(sim_bounded, package = "unitquantreg") sim_bounded_curr <- sim_bounded[sim_bounded$family == "uweibull", ] # Fitting using the analytical gradient fit_gradient <- unitquantreg(formula = y1 ~ x, data = sim_bounded_curr, tau = 0.5, family = "uweibull", control = unitquantreg.control(gradient = TRUE, trace = 1)) # Fitting without using the analytical gradient fit_nogradient <- unitquantreg(formula = y1 ~ x, data = sim_bounded_curr, tau = 0.5, family = "uweibull", control = unitquantreg.control(gradient = FALSE, trace = 1)) # Compare estimated coefficients cbind(gradient = coef(fit_gradient), no_gradient = coef(fit_nogradient))
data(sim_bounded, package = "unitquantreg") sim_bounded_curr <- sim_bounded[sim_bounded$family == "uweibull", ] # Fitting using the analytical gradient fit_gradient <- unitquantreg(formula = y1 ~ x, data = sim_bounded_curr, tau = 0.5, family = "uweibull", control = unitquantreg.control(gradient = TRUE, trace = 1)) # Fitting without using the analytical gradient fit_nogradient <- unitquantreg(formula = y1 ~ x, data = sim_bounded_curr, tau = 0.5, family = "uweibull", control = unitquantreg.control(gradient = FALSE, trace = 1)) # Compare estimated coefficients cbind(gradient = coef(fit_gradient), no_gradient = coef(fit_nogradient))
Density function, distribution function, quantile function and random number generation function
for the unit-Weibull distribution reparametrized in terms of the -th quantile,
.
duweibull(x, mu, theta, tau = 0.5, log = FALSE) puweibull(q, mu, theta, tau = 0.5, lower.tail = TRUE, log.p = FALSE) quweibull(p, mu, theta, tau = 0.5, lower.tail = TRUE, log.p = FALSE) ruweibull(n, mu, theta, tau = 0.5)
duweibull(x, mu, theta, tau = 0.5, log = FALSE) puweibull(q, mu, theta, tau = 0.5, lower.tail = TRUE, log.p = FALSE) quweibull(p, mu, theta, tau = 0.5, lower.tail = TRUE, log.p = FALSE) ruweibull(n, mu, theta, tau = 0.5)
x , q
|
vector of positive quantiles. |
mu |
location parameter indicating the |
theta |
nonnegative shape parameter. |
tau |
the parameter to specify which quantile use in the parametrization. |
log , log.p
|
logical; If TRUE, probabilities p are given as log(p). |
lower.tail |
logical; If TRUE, (default), |
p |
vector of probabilities. |
n |
number of observations. If |
Probability density function
Cumulative distribution function
Quantile function
Reparameterization
duweibull
gives the density, puweibull
gives the distribution function,
quweibull
gives the quantile function and ruweibull
generates random deviates.
Invalid arguments will return an error message.
Josmar Mazucheli
André F. B. Menezes
Mazucheli, J., Menezes, A. F. B and Ghitany, M. E., (2018). The unit-Weibull distribution and associated inference. Journal of Applied Probability and Statistics, 13(2), 1–22.
Mazucheli, J., Menezes, A. F. B., Fernandes, L. B., Oliveira, R. P. and Ghitany, M. E., (2020). The unit-Weibull distribution as an alternative to the Kumaraswamy distribution for the modeling of quantiles conditional on covariates. Journal of Applied Statistics, 47(6), 954–974.
Mazucheli, J., Menezes, A. F. B., Alqallaf, F. and Ghitany, M. E., (2021). Bias-Corrected Maximum Likelihood Estimators of the Parameters of the Unit-Weibull Distribution. Austrian Journal of Statistics, 50(3), 41–53.
set.seed(6969) x <- ruweibull(n = 1000, mu = 0.5, theta = 1.5, tau = 0.5) R <- range(x) S <- seq(from = R[1], to = R[2], by = 0.01) hist(x, prob = TRUE, main = 'unit-Weibull') lines(S, duweibull(x = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2) plot(ecdf(x)) lines(S, puweibull(q = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2) plot(quantile(x, probs = S), type = "l") lines(quweibull(p = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2)
set.seed(6969) x <- ruweibull(n = 1000, mu = 0.5, theta = 1.5, tau = 0.5) R <- range(x) S <- seq(from = R[1], to = R[2], by = 0.01) hist(x, prob = TRUE, main = 'unit-Weibull') lines(S, duweibull(x = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2) plot(ecdf(x)) lines(S, puweibull(q = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2) plot(quantile(x, probs = S), type = "l") lines(quweibull(p = S, mu = 0.5, theta = 1.5, tau = 0.5), col = 2)
Performs Vuong test between two fitted objects of class
unitquantreg
vuong.test(object1, object2, alternative = c("two.sided", "less", "greater"))
vuong.test(object1, object2, alternative = c("two.sided", "less", "greater"))
object1 , object2
|
objects of class |
alternative |
indicates the alternative hypothesis and must be one
of |
The statistic of Vuong likelihood ratio test for compare two non-nested regression models is defined by
where
is an estimator for the variance of
,
and
are the corresponding rival densities evaluated at the maximum likelihood estimates.
When we have that
in distribution.
Therefore, at
level of significance the null hypothesis of
the equivalence of the competing models is rejected if
,
where
is the
quantile of standard normal distribution.
In practical terms,
is better (worse) than
if
(or
).
A list with class "htest"
containing the following
components:
statistic |
the value of the test statistic. |
p.value |
the p-value of the test. |
alternative |
a character string describing the alternative hypothesis. |
method |
a character string with the method used. |
data.name |
a character string ginven the name of families models under comparison. |
André F. B. Menezes
Josmar Mazucheli
Vuong, Q. (1989). Likelihood ratio tests for model selection and non-nested hypotheses. Econometrica, 57(2), 307–333.
data(sim_bounded, package = "unitquantreg") sim_bounded_curr <- sim_bounded[sim_bounded$family == "uweibull", ] fit_uweibull <- unitquantreg(formula = y1 ~ x, tau = 0.5, data = sim_bounded_curr, family = "uweibull") fit_kum <- unitquantreg(formula = y1 ~ x, tau = 0.5, data = sim_bounded_curr, family = "kum") ans <- vuong.test(object1 = fit_uweibull, object2 = fit_kum) ans str(ans)
data(sim_bounded, package = "unitquantreg") sim_bounded_curr <- sim_bounded[sim_bounded$family == "uweibull", ] fit_uweibull <- unitquantreg(formula = y1 ~ x, tau = 0.5, data = sim_bounded_curr, family = "uweibull") fit_kum <- unitquantreg(formula = y1 ~ x, tau = 0.5, data = sim_bounded_curr, family = "kum") ans <- vuong.test(object1 = fit_uweibull, object2 = fit_kum) ans str(ans)
The access of people in households with piped water supply in the cities of Brazil from the Southeast and Northeast regions. Information obtained during the census of 2010.
data(water, package = "unitquantreg")
data(water, package = "unitquantreg")
data.frame
with 3457 observations and 5 columns:
phpws
: the proportion of households with piped water supply.
mhdi
: municipal human development index.
incpc
: per capita income.
region
: 0 for Southeast, 1 for Northeast.
pop
: total population.
André F. B. Menezes
Mazucheli, J., Menezes, A. F. B., Fernandes, L. B., Oliveira, R. P. and Ghitany, M. E., (2020). The unit-Weibull distribution as an alternative to the Kumaraswamy distribution for the modeling of quantiles conditional on covariates. Jounal of Applied Statistics, 47(6), 954–974.