Quantile-Based Estimator

Description

Quantile-based estimators (Q-estimators) can be used to fit any parametric distribution, using its quantile function. Q-estimators are usually more robust than standard maximum likelihood estimators. The method is described in: Sottile G. and Frumento P. (2022). Robust estimation and regression with parametric quantile functions. <doi:10.1016/j.csda.2022.107471>.

Details

The DESCRIPTION file:

Index of help topics:

Qcoxph                  Q-Estimation of Proportional Hazards Regression
                        Models
Qcoxph.control          Auxiliary for Controlling Qcoxph Fitting
Qest                    Q-Estimation
Qest-package            Quantile-Based Estimator
Qest.control            Auxiliary for Controlling Qest Fitting
Qfamily                 Family Objects for Qest
Qlm                     Q-Estimation of Linear Regression Models
Qlm.fit                 Fitter Functions for Quantile-based Linear
                        Models
invQ                    Inverse of Quantile Function
summary.Qest            Summarizing Q-estimators
wtrunc                  Weighting Function for 'Qest', 'Qlm', and
                        'Qcoxph'.

Author(s)

Gianluca Sottile [aut, cre], Paolo Frumento [aut]

Maintainer: Gianluca Sottile <gianluca.sottile@unipa.it>

References

Sottile G, and Frumento P (2022). Robust estimation and regression with parametric quantile functions. Computational Statistics and Data Analysis. <doi:10.1016/j.csda.2022.107471>

See Also

Qest, Qlm, Qcoxph

Examples

## Not run: 
Qest(y ~ x, Q, start) # General-purpose Q-estimator
Qlm(y ~ x) # Q-estimation of linear models
Qcoxph(Surv(time, event) ~ x) # Q-estimation of proportional hazards models

## End(Not run)

Q-Estimation of Proportional Hazards Regression Models

Description

Fit proportional hazards regression models using Q-estimation.

Usage

Qcoxph(formula, weights, start, data, knots, wtau = NULL,
    control = Qcoxph.control(), ...)

Arguments

Details

This function estimates a proportional hazards model, allowing for right-censored and left-truncated data. The syntax and output of Qcoxph are almost identical to those of coxph, but the parameters are estimated using Q-estimation (Sottile and Frumento, 2020). This method can be used to obtain outlier-robust estimators of the regression coefficients.

The quantile function of a proportional hazards model is given by

Q(\tau | x) = H0^{-1}(-exp{-x'\beta} log(1 - \tau))

where H0 is the baseline cumulative hazard function. In Qcoxph, H0 is parametrized by a piecewise linear function identified by the provided knots. As the number of knots increases, the baseline hazard becomes arbitrarily flexible.

Estimation is carried out by finding the zeroes of a set of integrals equation. The optional argument wtau permits assigning a different weight to each quantile in (0,1). It is possible to choose wtau to be a discontinuous function (e.g., wtau = function(tau){tau < 0.95}). However, this may occasionally result in poorly estimated of the standard errors.

The estimation algorithm is briefly described in the documentation of Qcoxph.control.

Value

an object of classes “Qcoxph”, “coxph”, and “Qest”. See coxph.object for details. All the S3 methods that are available for “coxph” objects will also work with a “Qcoxph” object.

An object of class “Qcoxph” is a list containing at least the following components:

Author(s)

Paolo Frumento <paolo.frumento@unipi.it>, Gianluca Sottile <gianluca.sottile@unipa.it>

References

Sottile G, and Frumento P (2022). Robust estimation and regression with parametric quantile functions. Computational Statistics and Data Analysis. <doi:10.1016/j.csda.2022.107471>

Muggeo VMR (2008). Segmented: an R package to fit regression models with broken-line relationships. R News 8/1, 20–25.

See Also

Qest, for general Q-estimation, and Qlm, for Q-estimation of linear models.

Examples

# A proportional-hazards Weibull model

n <- 100
x <- runif(n,0,3)
shape <- 2
t <- rweibull(n, shape = shape, scale = (1/exp(2 + 2*x))^(1/shape)) # time-to-event
c <- runif(n,0,1) # censoring variable
y <- pmin(t,c) # observed response
d <- (t <= c) # event indicator

require(survival)
m1 <- coxph(Surv(y,d) ~ x) # standard Cox model
m2 <- Qcoxph(Surv(y,d) ~ x) # Q-estimator

Auxiliary for Controlling Qcoxph Fitting

Description

Auxiliary function for controlling Qcoxph fitting. Estimation proceeds in three steps: (i) evaluation of starting points; (iia) stochastic gradient-based optimization (iib) standard gradient-based optimization; and (iii) Newton-Raphson. Step (i) is based on a preliminary fit of a Cox model via coxph. Steps (iia) and (iib) find an approximate solution, and make sure that the Jacobian matrix is well-defined. Finally, step (iii) finds a more precise solution.

Usage

Qcoxph.control(tol = 1e-8, maxit, safeit, alpha0, display = FALSE)

Arguments

Details

If called with no arguments, Qcoxph.control() returns a list with the current settings of these parameters. Any arguments included in the call sets those parameters to the new values, and then silently returns.

Value

A list with named elements as in the argument list

Author(s)

Gianluca Sottile <gianluca.sottile@unipa.it> Paolo Frumento <paolo.frumento@unipi.it>

See Also

Qcoxph

Q-Estimation

Description

An implementation of the quantile-based estimators described in Sottile and Frumento (2022).

Usage

Qest(formula, Q, weights, start, data, ntau = 199, wtau = NULL,
  control = Qest.control(), ...)

Arguments

Details

A parametric model, Q(\tau | \theta, x), is used to describe the conditional quantile function of an outcome Y, given a vector x of covariates. The model parameters, \theta, are estimated by minimizing the (weighted) integral, with respect to \tau, of the loss function of standard quantile regression. If the data are censored or truncated, \theta is estimated by solving a set of estimating equations. In either case, numerical integration is required to calculate the objective function: a grid of ntau points in (0,1) is used. The estimation algorithm is briefly described in the documentation of Qest.control.

The optional argument wtau can be used to attribute a different weight to each quantile. Although it is possible to choose wtau to be a discontinuous function (e.g., wtau = function(tau){tau < 0.95}), this may occasionally result in poorly estimated standard errors.

The quantile function Q must have at least the following three arguments: theta, tau, data, in this order. The first argument, theta, is a vector (not a matrix) of parameters' values. The second argument, tau, is the order of the quantile. When Q receives a n*ntau matrix of tau values, it must return a n*ntau matrix of quantiles. The third argument, data, is a data frame that includes the predictors used by Q.

If Q is identified by one Qfamily, everything becomes much simpler. It is not necessary to implement your own quantile function, and the starting points are not required. Note that ntau is ignored if Q = Qnorm or Q = Qunif.

Please check the documentation of Qfamily to see the available built-in distributions. A convenient Q-based implementation of the standard linear regression model is offered by Qlm. Proportional hazards models are implemented in Qcoxph.

Value

a list with the following elements:

Note

NOTE 1. If the data are censored or truncated, estimation is carried out by solving estimating equations, and no associated loss is present. In this case, a meaningful value of obj.function is the integrated loss [equation 1 of Sottile and Frumento (2022)] in which the indicator function I(y \le Q(\tau | \theta, x)) has been replaced with one of the expressions presented in equations 6 and 7 of the paper. The resulting loss, however, is not the function being minimized.

NOTE 2. To prevent computational problems, avoid situations in which some of the estimated parameters are expected to be very small or very large. For example, standardize the predictors, and normalize the response. Avoid as much as possible parameters with bounded support. For example, model a variance/rate/shape parameter on the log scale, e.g., \sigma = exp(\theta). Carefully select the starting points, and make sure that Q(start, ...) is well-defined. If Q is identified by one Qfamily, all these recommendations can be ignored.

NOTE 3. You should not use Qest to fit parametric models describing discrete distributions, where the quantile function is piecewise constant. You can try, but the optimization algorithm will most likely fail. The predefined family Qpois allows to fit a Poisson distribution by using a continuous version of its quantile function (see Qfamily).

Author(s)

Paolo Frumento <paolo.frumento@unipi.it>, Gianluca Sottile <gianluca.sottile@unipa.it>

References

Sottile G, and Frumento P (2022). Robust estimation and regression with parametric quantile functions. Computational Statistics and Data Analysis. <doi:10.1016/j.csda.2022.107471>

See Also

Qest.control, for operational parameters, and summary.Qest, for model summary. Qfamily, for the available built-in distributions. wtrunc for built-in weighting functions (wtau argument). Qlm, for Q-estimation of the standard normal (linear) regression model; Qcoxph, for proportional hazards models.

Examples

# Ex1. Normal model

# Quantile function of a linear model
Qlinmod <- function(theta, tau, data){
  sigma <- exp(theta[1])
  beta <- theta[-1]
  X <- model.matrix( ~ x1 + x2, data = data)
  qnorm(tau, X %*% beta, sigma)
}

n <- 100
x1 <- rnorm(n)
x2 <- runif(n,0,3)
theta <- c(1,4,1,2)
y <- Qlinmod(theta, runif(n), data.frame(x1,x2)) # generate the data

m1 <- Qest(y ~ x1 + x2, Q = Qlinmod, start = c(NA,NA,NA,NA)) # User-defined quantile function
summary(m1)

m2 <- Qest(y ~ x1 + x2, Q = Qnorm) # Qfamily
summary(m2)

m3 <- Qlm(y ~ x1 + x2)
summary(m3) # using 'Qlm' is much simpler and faster, with identical results





# Ex2. Weibull model with proportional hazards

# Quantile function
QWeibPH <- function(theta, tau, data){
  shape <- exp(theta[1])
  beta <- theta[-1]
  X <- model.matrix(~ x1 + x2, data = data)
  qweibull(tau, shape = shape, scale = (1/exp(X %*% beta))^(1/shape))
}

n <- 100
x1 <- rbinom(n,1,0.5)
x2 <- runif(n,0,3)
theta <- c(2,-0.5,1,1)

t <- QWeibPH(theta, runif(n), data.frame(x1,x2)) # time-to-event
c <- runif(n,0.5,1.5) # censoring variable
y <- pmin(t,c) # observed response
d <- (t <= c) # event indicator

m1 <- Qest(Surv(y,d) ~ x1 + x2, Q = QWeibPH, start = c(NA,NA,NA,NA))
summary(m1)

m2 <- Qcoxph(Surv(y,d) ~ x1 + x2)
summary(m2) # using 'Qcoxph' is much simpler and faster (but not identical)




# Ex3. A Gamma model

# Quantile function
Qgm <- function(theta, tau, data){
  a <- exp(theta[1])
  b <- exp(theta[2])
  qgamma(tau, shape = a, scale = b)
}
n <- 100
theta <- c(2,-1)
y <- rgamma(n, shape = exp(theta[1]), scale = exp(theta[2]))

m1 <- Qest(y ~ 1, Q = Qgm, start = c(NA, NA)) # User-defined quantile function
m2 <- Qest(y ~ 1, Q = Qgamma) # Qfamily
m3 <- Qest(y ~ 1, Q = Qgamma, wtau = function(tau, h) dnorm((tau - 0.5)/h), h = 0.2)
# In m3, more weight is assigned to quantiles around the median




# Ex4. A Poisson model

# Quantile function
n <- 100
x1 <- runif(n)
x2 <- rbinom(n,1,0.5)
y <- rpois(n, exp(1.5 -0.5*x1 + x2))
m1 <- Qest(y ~ x1 + x2, Q = Qpois) # Use a Qfamily! See "Note"
m2 <- Qest(y + runif(n) ~ x1 + x2, Q = Qpois) # Use jittering! See the documentation of "Qfamily"

Auxiliary for Controlling Qest Fitting

Description

Auxiliary function for controlling Qest fitting. Estimation proceeds in three steps: (i) evaluation of starting points; (iia) stochastic gradient-based optimization (iib) standard gradient-based optimization; and (iii) Newton-Raphson. Step (i) is initialized at the provided starting values (the start argument of Qest), and utilizes a preliminary flexible model, estimated with pchreg, to generate a cheap guess of the model parameters. If you have good starting points, you can skip step (i) by setting restart = FALSE. Steps (iia) and (iib) find an approximate solution, and make sure that the Jacobian matrix is well-defined. Finally, step (iii) finds a more precise solution.

Usage

Qest.control(tol = 1e-8, maxit, safeit, alpha0, display = FALSE, restart = FALSE)

Arguments

Details

If called with no arguments, Qest.control() returns a list with the current settings of these parameters. Any arguments included in the call sets those parameters to the new values, and then silently returns.

Value

A list with named elements as in the argument list

Note

Step (i) is not performed, and restart is ignored, if the quantile function is one of the available Qfamily.

Author(s)

Gianluca Sottile <gianluca.sottile@unipa.it> Paolo Frumento <paolo.frumento@unipi.it>

See Also

Qest and Qlm

Family Objects for Qest

Description

Family objects are used to specify the model to be fitted by Qest.

Usage

Qnorm()
Qgamma()
Qpois(offset = NULL)
Qunif(min = TRUE)

Arguments

Details

A Qfamily object can be used to identify a certain type of distribution within a call to Qest. You can supply either the name of the family, or the function itself, or a call to it. For example, the following are equivalent: Qest(formula, "Qpois"), Qest(formula, Qpois), and Qest(formula, Qpois()). The latter syntax can be used to pass additional arguments, if any.

The Qnorm family fits a normal homoskedastic model in which the mean is described by a linear predictor. The parameters are: log(sigma), beta. Qest(formula, Qnorm) is equivalent to Qlm(formula), but returns a very basic output. However, Qest allows for censored and truncated data, while Qlm does not.

The Qgamma family fits a Gamma distribution in which the log-scale is modeled by a linear predictor. The model parameters are: log(shape), beta.

The Qpois family fits a Poisson distribution in which the log-rate is modeled by a linear predictor. In reality, to obtain a continuous quantile function, qpois is replaced by the inverse, with respect to y, of the upper regularized gamma function, Q(y,\lambda). It is recommended to apply Qpois to a jittered response (i.e., y + runif(n)).

The Qunif family fits a Uniform distribution U(a,b) in which both a and b are modeled by linear predictors. The design matrix, however, is the same for a and b. Use Qunif(min = FALSE) to fit a U(0,b) model. The parameters are: beta_a, beta_b, or only beta_b if min = FALSE.

The families Qnorm and Qgamma can be used when the data are censored or truncated, while Qpois and Qunif cannot. All families can be estimated without covariates, using formula = ~ 1.

Value

An object of class "Qfamily" that contains all the necessary information to be passed to Qest.

Author(s)

Gianluca Sottile <gianluca.sottile@unipa.it>, Paolo Frumento <paolo.frumento@unipi.it>

See Also

Qest.

Examples

n <- 250
x <- runif(n)
eta <- 1 + 2*x # linear predictor

# Normal model
y <- rnorm(n, eta, exp(1))
m1 <- Qest(y ~ x, Qnorm)
# Use Qlm(y ~ x) instead!

# Gamma model
y <- rgamma(n, shape = exp(1), scale = exp(eta))
m2 <- Qest(y ~ x, Qgamma)

# Poisson model
y <- rpois(n, exp(eta))
m3 <- Qest(y ~ x, Qpois)
m4 <- Qest(y + runif(n) ~ x, Qpois) # Jittering is recommended

# Uniform model
y <- runif(n, 0, eta)
m5 <- Qest(y ~ x, Qunif(min = TRUE))  # U(a,b)
m6 <- Qest(y ~ x, Qunif(min = FALSE)) # U(0,b)

Q-Estimation of Linear Regression Models

Description

Use Q-estimation to fit a Normal model in which the mean is a linear function of the predictors, and the variance is constant.

Usage

Qlm(formula, data, subset, weights, na.action, start = NULL, contrasts = NULL,
  wtau = NULL, control = Qest.control(), ...)

Arguments

Details

This function is used exactly as lm, but estimates the model parameters as in Qest. Using Q-estimation allows to obtain outlier-robust estimators of the regression coefficients. The optional argument wtau permits assigning a different weight to each quantile in (0,1). It is possible to choose wtau to be a discontinuous function (e.g., wtau = function(tau){tau < 0.95}). However, this may occasionally result in poorly estimated of the standard errors.

Note that Qlm, like lm, does not support censored or truncated data.

Value

Qlm returns an object of classes “Qlm”, “lm”, and “Qest”. The generic accessor functions summary, coefficients, fitted.values, and residuals can be used to extract infromation from a “Qlm” object.

An object of class “Qlm” is a list containing at least the following elements:

Author(s)

Gianluca Sottile <gianluca.sottile@unipa.it>, Paolo Frumento <paolo.frumento@unipi.it>

References

Sottile G, and Frumento P (2022). Robust estimation and regression with parametric quantile functions. Computational Statistics and Data Analysis. <doi:10.1016/j.csda.2022.107471>

See Also

Qest, for general Q-estimation.

Examples

set.seed(1234)
n <- 100
x1 <- rnorm(n)
x2 <- runif(n,0,3)
theta <- c(1,4,1,2)
y <- rnorm(n, 4 + x1 + 2*x2, 1)

m1 <- Qlm(y ~ x1 + x2)
summary(m1)

Fitter Functions for Quantile-based Linear Models

Description

This is the basic computing engine called by “Qlm” used to fit quantile-based linear models. This function should only be used directly by experienced users.

Usage

Qlm.fit(y, X, w = rep(1, nobs), start = NULL, wtau = NULL,
  control = Qest.control(), ...)

Arguments

Value

a “list” with components

Author(s)

Gianluca Sottile <gianluca.sottile@unipa.it>, Paolo Frumento <paolo.frumento@unipi.it>

References

Sottile G, and Frumento P (2022). Robust estimation and regression with parametric quantile functions. Computational Statistics and Data Analysis. <doi:10.1016/j.csda.2022.107471>

See Also

Qlm

Examples

# Ex. 1 Normal model

set.seed(1234)
n <- 100
x1 <- rnorm(n)
x2 <- runif(n,0,3)
y <- rnorm(n, 4 + x1 + 2*x2, 1)
X <- cbind(1, x1, x2)
w <- rep.int(1, n)

m <- Qlm.fit(y = y, X = X, w = w, control = Qest.control(display = TRUE))

Internal Functions

Description

Functions for internal use only, or not yet documented.

Usage

expit(x)
logit(x)
pmax0(x)
num.fun(dx,fx)
formatPerc(probs, digits)
Ltau(opt, tau)
minabs(x1,x2)
invJ(J, type)
tensorX(X)
buildTau(ntau, wtau = NULL, nobs, wtauoptions = NULL)
callwtau(wtau, tau, opt)
callQ(Q, theta, tau, data)
start.Qest(z, y, d, x, w, tau, Q, opt, start, data, type, control)
start.Qest.family(z, y, d, x, w, tau, wtau, wtauoptions,
  Q, opt, start, data, type, control)
rq.fit.br2(x, y, tau = 0.5)

Qest.sgs.internal(theta, type, tol, maxit, alpha0, ee, display, eps, n.it, ...)
Qest.gs.internal(theta, type, tol, maxit, alpha0, ee, display, eps, n.it, ...)
Qest.gs(theta, type, tol, maxit, alpha0, ee, display, eps, ...)
Qest.newton(theta, type, tol, maxit, safeit, alpha0, display, eps, ...)

Qlm.bfun(wtau, ...)
start.Qlm(x, y, w, start, ok, Stats)
scalevars.Qlm(X,y)
descalecoef.Qlm(theta, Stats)
Qlm.sgs.internal(theta, type, tol, maxit, alpha0, ee, display, n.it, y, X, w, bfun)
Qlm.gs.internal(theta, type, tol, maxit, alpha0, ee, display, n.it, y, X, w, bfun)
Qlm.gs(theta, type, tol, maxit, alpha0, ee, display, y, X, w, bfun)
Qlm.newton(theta, type = "u", tol, maxit, safeit, alpha0, display, y, X, w, bfun)

plfcox(y, knots, deriv = 0)
scalevars.Qcoxph(X,z,y,knots)
descalecoef.Qcoxph(theta, Stats)
check.singularities(X, scaleVars)
starting.points.Qcox(X, Y, n, w, mf, knots)
adjust.coef(theta)
agsurv.Qcoxph(y, x, wt, risk, fit)
basehaz.Qcoxph(fit, centered = TRUE, se.fit = FALSE)
coxsurv.fit.Qcoxph(ctype, stype, se.fit, varmat, cluster,
  y, x, wt, risk, position, strata, oldid, y2, x2, risk2,
  strata2, id2, unlist = TRUE, fit)
seg.lm.fit1(y,XREG,Z,PSI,return.all.sol=FALSE)

gs(theta0, f, ..., tol = 1e-4, maxit = 100)
myg(theta, f, f0, eps, ...)
dlist(x1,x2)
omega(d, tau, type, Fy, Fz)
choose_eps(Q, theta, y, data, eps0, obj = 0.01)
derQtheta(theta, eps, Q, Q1, data, tau, ind)
der2Qtheta(theta, eps, Q, Qtheta1, data, tau)
intA(A, tau, ifun = TRUE)
findp(y, tau, Q1)
findAp(p,tau,A)

A_beta_fun(BB)
A_gamma_fun(BB)
A_beta_beta_fun(BB)
A_gamma_gamma_mix_fun(BB)
A_gamma_gamma_fun(BB)
A_beta_gamma_fun(BB)
coxBB(theta, y, X, knots, tau)

derQtheta.gamma(Q)
der2Qtheta.gamma(Q, Qtheta)
findAp.gamma(atau, tau, dtau, p, int = TRUE)
QestGamma.ee.u(theta, eps, z, y, d, Q, w, data, tau, J = FALSE, EE)
QestGamma.ee.c(theta, eps, z, y, d, Q, w, data, tau, J = FALSE, EE)
QestGamma.ee.ct(theta, eps, z, y, d, Q, w, data, tau, J = FALSE, EE)

tau.pois(tau)
ppoisC(y, lambda)
dpoisC(y, lambda)
qpoisC.955(z, lambda)
qpoisC.me(log.lambda, A, B)
qpoisC.bisec(tau, lambda)
qpoisC(obj)

derQtheta.pois(Q)
der2Qtheta.pois(Q, Qtheta)
findp.pois(y, tau, Q1, Fy, theta)
findAp.pois(p, tau, A)
QestPois.ee.u(theta, eps, z, y, d, Q, w, data, tau, J = FALSE, EE)

QestUnif.ee.u(theta, eps, z, y, d, Q, w, data, tau, J = FALSE, EE)

QestNorm.ee.u(theta, eps, z, y, d, Q, w, data, tau, J = FALSE, EE)
QestNorm.ee.c(theta, eps, z, y, d, Q, w, data, tau, J = FALSE, EE)
QestNorm.ee.ct(theta, eps, z, y, d, Q, w, data, tau, J = FALSE, EE)

Qest.ee.u(theta, eps, z, y, d, Q, w, data, tau, J = FALSE, EE)
Qest.ee.c(theta, eps, z, y, d, Q, w, data, tau, J = FALSE, EE)
Qest.ee.ct(theta, eps, z, y, d, Q, w, data, tau, J = FALSE, EE)
QCox.ee.c(theta, eps, z, y, d, X, w, knots, tau, J = FALSE, EE)
QCox.ee.ct(theta, eps, z, y, d, X, w, knots, tau, J = FALSE, EE)
Qlm.ee.u(theta, X, w, bfun, EE, J = FALSE)

Qest.covar(fit, eps, w)
Qcox.covar(theta, z, y, d, X, w, knots, tau, type)
Qlm.covar(g.i, w, H)

Loss(w, d, tau, type, Fy, Fz)
coxLoss(theta, z, y, d, X, w, knots, tau, type, Fy, Fz)
qlmLoss(theta, y, X, w, bfun)

## S3 method for class 'Qest'
print(x, digits = max(3L, getOption("digits") - 3L), ...)
## S3 method for class 'summary.Qest'
print(x, digits = max(3L, getOption("digits") - 3L), ...)
## S3 method for class 'Qest'
confint(object, parm, level = 0.95, ...)
## S3 method for class 'Qest'
vcov(object, ...)

## S3 method for class 'Qlm'
summary(object, correlation = FALSE,
  symbolic.cor = FALSE, ...)
## S3 method for class 'summary.Qlm'
print(x, digits = max(3L, getOption("digits") - 3L),
  symbolic.cor = x$symbolic.cor, signif.stars = getOption("show.signif.stars"),
  ...)
## S3 method for class 'Qlm'
vcov(object, ...)

## S3 method for class 'Qcoxph'
print(x, digits = max(1L, getOption("digits") - 3L),
  signif.stars = FALSE, ...)
## S3 method for class 'Qcoxph'
summary(object, conf.int = 0.95, scale = 1, ...)
## S3 method for class 'summary.Qcoxph'
print(x, digits = max(getOption("digits") - 3, 3),
  signif.stars = getOption  ("show.signif.stars"), ...)
## S3 method for class 'Qcoxph'
survfit(formula, newdata, se.fit = TRUE, conf.int = 0.95,
  individual = FALSE, stype = 2, ctype, conf.type = c("log", "log-log",
  "plain","none", "logit", "arcsin"), censor = TRUE, start.time, id,
  influence = FALSE, na.action = na.pass, type, ...)
## S3 method for class 'Qcoxph'
residuals(object, type = c("martingale", "deviance", "score",
  "schoenfeld", "dfbeta", "dfbetas", "scaledsch", "partial"),
  collapse = FALSE, weighted = FALSE, ...)
## S3 method for class 'Qcoxph'
predict(object, newdata, type = c("lp", "risk", "expected",
  "terms", "survival"), se.fit = FALSE, na.action = na.pass,
  terms = names(object$assign), collapse, reference = c("strata", "sample"),
  ...)

Value

No return value, internal functions.

Inverse of Quantile Function

Description

Auxiliary function to compute cumulative distribution function (CDF) by inverting a quantile function.

Usage

invQ(Q, theta, y, data, n.it = 17)

Arguments

Details

Given a parametric quantile function Q(\tau | \theta), the CDF is defined as F(y | \theta) = Q^{-1}(y | \theta). Alternatively, F(y | \theta) corresponds to the value \tau* such that Q(\tau* | \theta) = y. Starting from \tau = 0.5, a bisection algorithm is used to evaluate numerically \tau*. The maximum error is given by 1/2^(n.it + 1).

Value

a vector of CDF values.

Author(s)

Maintainer: Gianluca Sottile <gianluca.sottile@unipa.it>

See Also

Qest

Examples

# Ex. 1 Normal model

# Quantile function of a linear model.
Qlinmod <- function(theta, tau, data){
  sigma <- exp(theta[1])
  beta <- theta[-1]
  X <- model.matrix( ~ x1 + x2, data = data)
  qnorm(tau, X %*% beta, sigma)
}

n <- 100
x1 <- rnorm(n)
x2 <- runif(n,0,3)
theta <- c(1,4,1,2)

# generate the data
U <- runif(n)
y <- Qlinmod(theta, U, data.frame(x1,x2))

# Given y and theta, evaluate U = F(y)
invQ(Qlinmod, theta, y, data.frame(x1,x2))

Summarizing Q-estimators

Description

Summary method for class “Qest”.

Usage

## S3 method for class 'Qest'
summary(object, covar = FALSE, ...)

Arguments

Details

This function returns a summary of the most relevant information on model parameters, standard errors, and convergence status.

Value

The function summary.Qest computes and returns a list of summary statistics of the fitted model given in object, using the "call" and "terms" from its argument, plus

Author(s)

Gianluca Sottile <gianluca.sottile@unipa.it>

References

Sottile G, and Frumento P (2022). Robust estimation and regression with parametric quantile functions. Computational Statistics and Data Analysis. <doi:10.1016/j.csda.2022.107471>

See Also

Qest, for model fitting.

Examples

# Quantile function of an Exponential model
Qexp <- function(theta, tau, data){
  qexp(tau, exp(theta))
}

y <- rexp(100, exp(1))
m1 <- Qest(y ~ 1, Q = Qexp, start = NA)
summary(m1)
summary(m1, covar = TRUE)

Weighting Function for Qest, Qlm, and Qcoxph.

Description

This function can be used within a call to Qest, Qlm, or Qcoxph to assign a different weight to each quantile.

Usage

wtrunc(tau, delta.left = 0.05, delta.right = 0.05, smooth = FALSE, sigma = 0.01)

Arguments

Details

Within a call to Qest, Qlm, or Qcoxph, one may want to assign a different weight to each quantile through the optional argument wtau. This can be done for reasons of efficiency, or to counteract the presence of outliers. While wtau can be any user-defined function, one can use wtrunc as a shortcut to construct a weighting function that truncates a certain subset of quantiles in the tails of the distribution. For instance, the estimator defined by Qest(..., wtau = wtrunc, delta.left = 0.05, delta.right = 0.1) only uses quantiles in the interval (0.05, 0.90) to fit the model. In this example, delta.left = 0.05 is a guess for the proportion of outliers on the left tail; and delta.right is a guess for the proportion of outliers on the right tail. Use smooth = TRUE to replace the indicator functions involved in wtrunc with smooth functions. Introducing a weighting function that only assigns a positive weight to the quantiles that correspond to the “healthy” part of the distribution allows to deal with any level of contamination by outliers.

Value

A vector of weights assigned to each quantile.

Author(s)

Gianluca Sottile <gianluca.sottile@unipa.it>, Paolo Frumento <paolo.frumento@unipi.it>

See Also

Qest, Qlm, Qcoxph.

Examples

## Not run: 
taus <- seq(0, 1, length.out = 1000)

### zero weight to quantiles above 0.95
plot(taus, wtrunc(taus, delta.left = 0, delta.right = 0.05),
  type = "l", lwd = 1.5)
# smooth counterpart
lines(taus, wtrunc(taus, delta.left = 0, delta.right = 0.05,
  smooth = TRUE, sigma = .01), col = 2, lwd = 1.5)
lines(taus, wtrunc(taus, delta.left = 0, delta.right = 0.05,
  smooth = TRUE, sigma = .05), col = 3, lwd = 1.5)

### zero weight to quantiles below 0.05
plot(taus, wtrunc(taus, delta.left = 0.05, delta.right = 0),
  type = "l", lwd = 1.5)
# smooth counterpart
lines(taus, wtrunc(taus, delta.left = 0.05, delta.right = 0,
  smooth = TRUE, sigma = .01), col = 2, lwd = 1.5)
lines(taus, wtrunc(taus, delta.left = 0.05, delta.right = 0,
  smooth = TRUE, sigma = .05), col = 3, lwd = 1.5)


### zero weight to quantiles below 0.05 and above 0.90
plot(taus, wtrunc(taus, delta.left = 0.05, delta.right = 0.10),
  type = "l", lwd = 1.5)
# smooth counterpart
lines(taus, wtrunc(taus, delta.left = 0.05, delta.right = 0.10,
  smooth = TRUE, sigma = .01), col = 2, lwd = 1.5)
lines(taus, wtrunc(taus, delta.left = 0.05, delta.right = 0.10,
  smooth = TRUE, sigma = .05), col = 3, lwd = 1.5)


### Use wtrunc in Qest, Qlm, Qcoxph
# Qest(..., wtau = wtrunc, delta.left = 0.05, delta.right = 0.1)

## End(Not run)