amlnormal: Asymmetric Least Squares Quantile Regression

Description

Asymmetric least squares, a special case of maximizing an asymmetric likelihood function of a normal distribution. This allows for expectile/quantile regression using asymmetric least squares error loss.

Usage

amlnormal(w.aml = 1, parallel = FALSE, lexpectile = "identity",
          eexpectile = list(), iexpectile = NULL, imethod = 1, digw = 4)

Arguments

w.aml

Numeric, a vector of positive constants controlling the percentiles. The larger the value the larger the fitted percentile value (the proportion of points below the ``w-regression plane''). The default value of unity results in the ordinary least sq

parallel

If w.aml has more than one value then this argument allows the quantile curves to differ by the same amount as a function of the covariates. Setting this to be TRUE should force the quantile curves to not cross (although

lexpectile, eexpectile, iexpectile

See CommonVGAMffArguments for more information.

imethod

Integer, either 1 or 2 or 3. Initialization method. Choose another value if convergence fails.

digw

Passed into Round as the digits argument for the w.aml values; used cosmetically for labelling.

Value

An object of class "vglmff" (see vglmff-class). The object is used by modelling functions such as vglm and vgam.

Details

This is an implementation of Efron (1991) and full details can be obtained there. Equation numbers below refer to that article. The model is essentially a linear model (see lm), however, the asymmetric squared error loss function for a residual $r$ is $r^2$ if $r \leq 0$ and $w r^2$ if $r > 0$. The solution is the set of regression coefficients that minimize the sum of these over the data set, weighted by the weights argument (so that it can contain frequencies). Newton-Raphson estimation is used here.

References

Efron, B. (1991) Regression percentiles using asymmetric squared error loss. Statistica Sinica, 1, 93--125.

Examples

Run this code

# Example 1
ooo = with(bminz, order(age))
bminz = bminz[ooo,]  # Sort by age
(fit = vglm(BMI ~ bs(age), fam=amlnormal(w.aml=0.1), bminz))
fit@extra  # Gives the w value and the percentile
coef(fit, matrix=TRUE)

# Quantile plot
with(bminz, plot(age, BMI, col="blue", main=
     paste(round(fit@extra$percentile, dig=1),
           "expectile-percentile curve")))
with(bminz, lines(age, c(fitted(fit)), col="black"))



# Example 2
# Find the w values that give the 25, 50 and 75 percentiles
findw = function(w, percentile=50) {
    fit2 = vglm(BMI ~ bs(age), fam=amlnormal(w=w), data=bminz)
    fit2@extra$percentile - percentile
}
# Quantile plot
with(bminz, plot(age, BMI, col="blue", las=1, main=
     "25, 50 and 75 expectile-percentile curves"))
for(myp in c(25,50,75)) {
# Note: uniroot() can only find one root at a time
    bestw = uniroot(f=findw, interval=c(1/10^4, 10^4), percentile=myp)
    fit2 = vglm(BMI ~ bs(age), fam=amlnormal(w=bestw$root), data=bminz)
with(bminz, lines(age, c(fitted(fit2)), col="red"))
}



# Example 3; this is Example 1 but with smoothing splines and
# a vector w and a parallelism assumption.
ooo = with(bminz, order(age))
bminz = bminz[ooo,]  # Sort by age
fit3 = vgam(BMI ~ s(age, df=4), fam=amlnormal(w=c(.1,1,10), parallel=TRUE),
            bminz, trac=TRUE)
fit3@extra # The w values, percentiles and weighted deviances

# The linear components of the fit; not for human consumption:
coef(fit3, matrix=TRUE)

# Quantile plot
with(bminz, plot(age, BMI, col="blue", main=
     paste(paste(round(fit3@extra$percentile, dig=1), collapse=", "),
           "expectile-percentile curves")))
with(bminz, matlines(age, fitted(fit3), col=1:fit3@extra$M, lwd=2))
with(bminz, lines(age, c(fitted(fit )), col="black")) # For comparison

Run the code above in your browser using DataLab

The Learning Leader's Guide to AI Literacy