alsqreg {VGAM} | R Documentation |
Quantile regression using asymmetric least squares error loss.
alsqreg(w.als=1, parallel=FALSE, lexpectile = "identity", eexpectile = list(), iexpectile = NULL, method.init=1, digw=4)
w.als |
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.als 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, eexpectile, iexpectile |
See CommonVGAMffArguments for more information.
|
method.init |
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.als values;
used cosmetically for labelling.
|
This method was proposed by 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.
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.als"
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.als
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 alsqreg
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.als
has more than one value then the value returned by
the slot is the sum taken over all the w.als
values.
This VGAM family function could well be renamed amlnormal()
instead, given the other function names amlpoisson
,
amlbinomial
, etc.
Thomas W. Yee
Efron, B. (1991) Regression percentiles using asymmetric squared error loss. Statistica Sinica, 1, 93–125.
amlpoisson
,
amlbinomial
,
amlexponential
,
bminz
,
alaplace1
,
lms.bcn
and similar variants are alternative
methods for quantile regression.
# Example 1 data(bminz) o = with(bminz, order(age)) bminz = bminz[o,] # Sort by age (fit = vglm(BMI ~ bs(age), fam=alsqreg(w.als=0.1), data=bminz)) fit@extra # Gives the w value and the percentile coef(fit) coef(fit, matrix=TRUE) ## Not run: # 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")) ## End(Not run) # 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=alsqreg(w=w), data=bminz) fit2@extra$percentile - percentile } ## Not run: # Quantile plot with(bminz, plot(age, BMI, col="blue", las=1, main= "25, 50 and 75 expectile-percentile curves")) ## End(Not run) 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=alsqreg(w=bestw$root), data=bminz) ## Not run: with(bminz, lines(age, c(fitted(fit2)), col="red")) ## End(Not run) } # Example 3; this is Example 1 but with smoothing splines and # a vector w and a parallelism assumption. data(bminz) o = with(bminz, order(age)) bminz = bminz[o,] # Sort by age fit3 = vgam(BMI ~ s(age, df=4), fam=alsqreg(w=c(.1,1,10), parallel=TRUE), data=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) ## Not run: # 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 ## End(Not run)