# amlnormal: Asymmetric Least Squares Quantile Regression In VGAM: Vector Generalized Linear and Additive Models

## 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).

Thomas W. Yee

## References

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

`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

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.