amlnormal | R Documentation |
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.
amlnormal(w.aml = 1, parallel = FALSE, lexpectile = "identitylink",
iexpectile = NULL, imethod = 1, digw = 4)
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 |
lexpectile , iexpectile |
See |
imethod |
Integer, either 1 or 2 or 3. Initialization method. Choose another value if convergence fails. |
digw |
Passed into |
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.
An object of class "vglmff"
(see
vglmff-class
). The object is used by modelling
functions such as vglm
and vgam
.
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
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.
## 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), 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), 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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.