amlnormal: Asymmetric Least Squares Quantile Regression

Description Usage Arguments Details Value Note Author(s) References See Also Examples

View source: R/family.qreg.R

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

1
2
amlnormal(w.aml = 1, parallel = FALSE, lexpectile = "identitylink",
          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 squares (OLS) solution.

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 they may not cross anyway). See CommonVGAMffArguments for more information.

lexpectile, 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.

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 <= 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.

Value

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

Note

On fitting, the extra slot has list components "w.aml" and "percentile". The latter is the percent of observations below the “w-regression plane”, which is the fitted values.

One difficulty is finding the w.aml value giving a specified percentile. One solution is to fit the model within a root finding function such as uniroot; see the example below.

For amlnormal objects, methods functions for the generic functions qtplot and cdf have not been written yet.

See the note in amlpoisson on the jargon, including expectiles and regression quantiles.

The deviance slot computes the total asymmetric squared error loss (2.5). If w.aml has more than one value then the value returned by the slot is the sum taken over all the w.aml values.

This VGAM family function could well be renamed amlnormal() instead, given the other function names amlpoisson, amlbinomial, etc.

In this documentation the word quantile can often be interchangeably replaced by expectile (things are informal here).

Author(s)

Thomas W. Yee

References

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

See Also

amlpoisson, amlbinomial, amlexponential, bmi.nz, extlogF1, alaplace1, denorm, lms.bcn and similar variants are alternative methods for quantile regression.

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
## Not run: 
# Example 1
ooo <- with(bmi.nz, order(age))
bmi.nz <- bmi.nz[ooo, ]  # Sort by age
(fit <- vglm(BMI ~ sm.bs(age), amlnormal(w.aml = 0.1), data = bmi.nz))
fit@extra  # Gives the w value and the percentile
coef(fit, matrix = TRUE)

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

# Example 2
# Find the w values that give the 25, 50 and 75 percentiles
find.w <- function(w, percentile = 50) {
  fit2 <- vglm(BMI ~ sm.bs(age), amlnormal(w = w), data = bmi.nz)
  fit2@extra$percentile - percentile
}
# Quantile plot
with(bmi.nz, 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 = find.w, interval = c(1/10^4, 10^4), percentile = myp)
  fit2 <- vglm(BMI ~ sm.bs(age), amlnormal(w = bestw$root), data = bmi.nz)
  with(bmi.nz, lines(age, c(fitted(fit2)), col = "orange"))
}

# Example 3; this is Example 1 but with smoothing splines and
# a vector w and a parallelism assumption.
ooo <- with(bmi.nz, order(age))
bmi.nz <- bmi.nz[ooo, ]  # Sort by age
fit3 <- vgam(BMI ~ s(age, df = 4), data = bmi.nz, trace = TRUE,
             amlnormal(w = c(0.1, 1, 10), parallel = 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(bmi.nz, plot(age, BMI, col="blue", main =
     paste(paste(round(fit3@extra$percentile, digits = 1), collapse = ", "),
           "expectile-percentile curves")))
with(bmi.nz, matlines(age, fitted(fit3), col = 1:fit3@extra$M, lwd = 2))
with(bmi.nz, lines(age, c(fitted(fit )), col = "black"))  # For comparison

## End(Not run)

Example output

Loading required package: stats4
Loading required package: splines

Call:
vglm(formula = BMI ~ sm.bs(age), family = amlnormal(w.aml = 0.1), 
    data = bmi.nz)


Coefficients:
(Intercept) sm.bs(age)1 sm.bs(age)2 sm.bs(age)3 
  22.492304   -1.114858    7.496455   -5.124937 

Degrees of Freedom: 700 Total; 696 Residual
Residual deviance: 2986.426 
$w.aml
[1] 0.1

$M
[1] 1

$n
[1] 700

$y.names
[1] "w.aml = 0.1"

$percentile
w.aml = 0.1 
   20.14286 

$individual
[1] TRUE

$deviance
w.aml = 0.1 
   2986.426 

            expectile(w.aml = 0.1)
(Intercept)              22.492304
sm.bs(age)1              -1.114858
sm.bs(age)2               7.496455
sm.bs(age)3              -5.124937
Warning messages:
1: In vglm.fitter(x = x, y = y, w = w, offset = offset, Xm2 = Xm2,  :
  iterations terminated because half-step sizes are very small
2: In vglm.fitter(x = x, y = y, w = w, offset = offset, Xm2 = Xm2,  :
  some quantities such as z, residuals, SEs may be inaccurate due to convergence at a half-step
VGAM  s.vam  loop  1 :  deviance = 77690.291
VGAM  s.vam  loop  2 :  deviance = 72019.671
VGAM  s.vam  loop  3 :  deviance = 71819.18
VGAM  s.vam  loop  4 :  deviance = 71807.834
VGAM  s.vam  loop  5 :  deviance = 71807.656
VGAM  s.vam  loop  6 :  deviance = 71807.656
$w.aml
[1]  0.1  1.0 10.0

$M
[1] 3

$n
[1] 700

$y.names
[1] "w.aml = 0.1" "w.aml = 1"   "w.aml = 10" 

$percentile
w.aml = 0.1   w.aml = 1  w.aml = 10 
   19.85714    56.85714    85.57143 

$individual
[1] TRUE

$deviance
w.aml = 0.1   w.aml = 1  w.aml = 10 
   2972.058   14284.349   54551.249 

               expectile(w.aml = 0.1) expectile(w.aml = 1)
(Intercept)               22.76123490          26.15205226
s(age, df = 4)             0.01246891           0.01246891
               expectile(w.aml = 10)
(Intercept)              30.54240967
s(age, df = 4)            0.01246891

VGAM documentation built on Jan. 16, 2021, 5:21 p.m.