alaplace3: Asymmetric Laplace Distribution Family Functions
In VGAMdata: Data Supporting the 'VGAM' Package
View source: R/family.laplace.R
alaplace R Documentation
Asymmetric Laplace Distribution Family Functions
Description
Maximum likelihood estimation of
the 1, 2 and 3-parameter asymmetric Laplace distributions (ALDs).
The 2-parameter ALD may,
with trepidation and lots of skill,
sometimes be used as an approximation of quantile regression.
Usage
alaplace1(tau = NULL, llocation = "identitylink",
ilocation = NULL, kappa = sqrt(tau/(1 - tau)), Scale.arg = 1,
ishrinkage = 0.95, parallel.locat = TRUE ~ 0, digt = 4,
idf.mu = 3, zero = NULL, imethod = 1)
alaplace2(tau = NULL, llocation = "identitylink", lscale = "loglink",
ilocation = NULL, iscale = NULL, kappa = sqrt(tau/(1 - tau)),
ishrinkage = 0.95,
parallel.locat = TRUE ~ 0,
parallel.scale = FALSE ~ 0,
digt = 4, idf.mu = 3, imethod = 1, zero = "scale")
alaplace3(llocation = "identitylink", lscale = "loglink",
lkappa = "loglink", ilocation = NULL, iscale = NULL,
ikappa = 1, imethod = 1, zero = c("scale", "kappa"))
Arguments
tau, kappa
Numeric vectors with
0 < \tau < 1 and \kappa >0.
Most users will only specify tau since the estimated
location parameter corresponds to the \tauth
regression quantile, which is easier to understand.
See below for details.
llocation, lscale, lkappa
Character.
Parameter link functions for
location parameter \xi,
scale parameter \sigma,
asymmetry parameter \kappa.
See Links for more choices.
For example, the argument llocation can help handle
count data by restricting the quantiles to be positive
(use llocation = "loglink").
However, llocation is best left alone since the theory
only works properly with the identity link.
ilocation, iscale, ikappa
Optional initial values.
If given, it must be numeric and values are recycled to the
appropriate length.
The default is to choose the value internally.
parallel.locat, parallel.scale
See the parallel argument
of CommonVGAMffArguments.
These arguments apply to the location and scale parameters.
It generally only makes sense for the scale parameters
to be equal, hence set parallel.scale = TRUE.
Note that
assigning parallel.locat the value TRUE
circumvents the seriously embarrassing quantile crossing
problem because all constraint matrices except for the intercept
correspond to a parallelism assumption.
imethod
Initialization method.
Either the value 1, 2, 3 or 4.
idf.mu
Degrees of freedom for the cubic smoothing spline fit
applied to get an initial estimate of the location parameter.
See vsmooth.spline.
Used only when imethod = 3.
ishrinkage
How much shrinkage is used when initializing \xi.
The value must be between 0 and 1 inclusive, and
a value of 0 means the individual response values are used,
and a value of 1 means the median or mean is used.
This argument is used only when imethod = 4.
See CommonVGAMffArguments for more information.
Scale.arg
The value of the scale parameter \sigma.
This argument may be used to compute quantiles at
different \tau values from an existing fitted
alaplace2() model (practical only if it has a
single value).
If the model has parallel.locat = TRUE then
only the intercept need be estimated; use an offset.
See below for an example.
digt
Passed into Round as the digits
argument for the tau values; used cosmetically for
labelling.
zero
See CommonVGAMffArguments for more information.
Where possible,
the default is to model all the \sigma
and \kappa as an intercept-only term.
See CommonVGAMffArguments for more information.
Details
These VGAM family functions implement one variant of asymmetric
Laplace distributions (ALDs) suitable for quantile regression.
Kotz et al. (2001) call it the ALD.
Its density function is
f(y;\xi,\sigma,\kappa) = \frac{\sqrt{2}}{\sigma} \,
\frac{\kappa}{1 + \kappa^2} \,
\exp \left( - \frac{\sqrt{2}}{\sigma \, \kappa} |y - \xi |
\right)
for y \leq \xi, and
f(y;\xi,\sigma,\kappa) = \frac{\sqrt{2}}{\sigma} \,
\frac{\kappa}{1 + \kappa^2} \,
\exp \left( - \frac{\sqrt{2} \, \kappa}{\sigma} |y - \xi |
\right)
for y > \xi.
Here, the ranges are for all real y and
\xi, positive \sigma
and positive \kappa.
The special case \kappa = 1 corresponds to the
(symmetric) Laplace distribution of Kotz et al. (2001).
The mean is \xi + \sigma (1/\kappa - \kappa) / \sqrt{2}
and the variance is
\sigma^2 (1 + \kappa^4) / (2 \kappa^2).
The enumeration of the linear/additive predictors used for
alaplace2() is
the first location parameter followed by the first scale parameter,
then the second location parameter followed by the
second scale parameter, etc.
For alaplace3(), only a vector response is handled
and the last (third) linear/additive predictor is for
the asymmetry parameter.
It is known that the maximum likelihood estimate of the
location parameter \xi corresponds to the regression
quantile estimate of the classical quantile regression approach
of Koenker and Bassett (1978). An important property of the
ALD is that
P(Y \leq \xi) = \tau where
\tau = \kappa^2 / (1 + \kappa^2)
so that
\kappa = \sqrt{\tau / (1-\tau)}.
Thus alaplace2() might be used as an alternative to rq
in the quantreg package, although scoring is really
an unsuitable algorithm for estimation here.
Both alaplace1() and alaplace2() can handle
multiple responses, and the number of linear/additive
predictors is dictated by the length of tau or
kappa. The functions alaplace1()
and alaplace2() can also
handle multiple responses (i.e., a matrix response)
but only with a single-valued tau or kappa.
Value
An object of class "vglmff" (see
vglmff-class). The object is used by modelling
functions such as vglm and vgam.
In the extra slot of the fitted object are some list
components which are useful, e.g., the sample proportion of
values which are less than the fitted quantile curves.
Warning
These functions are experimental and especially subject to
change or withdrawal.
The usual MLE regularity conditions do not hold for
this distribution so that misleading inferences may result,
e.g., in the summary and vcov of the object. The
1-parameter ALD can be approximated by extlogF1
which has continuous derivatives and is recommended over
alaplace1.
Care is needed with tau values which are too small,
e.g., for count data with llocation = "loglink" and if
the sample proportion of zeros is greater than tau.
Note
These VGAM family functions use Fisher scoring.
Convergence may be slow and half-stepping is usual
(although one can use trace = TRUE to see which
is the best model and then use maxit to choose
that model)
due to the regularity conditions not holding.
Often the iterations slowly crawl towards the solution so
monitoring the convergence (set trace = TRUE) is highly
recommended.
Instead, extlogF1 is recommended.
For large data sets it is a very good idea to keep the length of
tau/kappa low to avoid large memory requirements.
Then
for parallel.locat = FALSE one can repeatedly fit a model with
alaplace1() with one \tau at a time;
and
for parallel.locat = TRUE one can refit a model with
alaplace1() with one \tau at a time but
using offsets and an intercept-only model.
A second method for solving the noncrossing quantile problem is
illustrated below in Example 3.
This is called the accumulative quantile method (AQM)
and details are in Yee (2015).
It does not make the strong parallelism assumption.
The functions alaplace2() and laplace
differ slightly in terms of the parameterizations.
Author(s)
Thomas W. Yee
References
Koenker, R. and Bassett, G. (1978).
Regression quantiles.
Econometrica,
46, 33–50.
Kotz, S., Kozubowski, T. J. and Podgorski, K. (2001).
The Laplace distribution and generalizations:
a revisit with applications to communications,
economics, engineering, and finance,
Boston: Birkhauser.
See Also
ralap,
laplace,
extlogF1,
CommonVGAMffArguments,
lms.bcn,
amlnormal,
sc.studentt2,
simulate.vlm.
Examples
## Not run:
# Example 1: quantile regression with smoothing splines
set.seed(123); adata <- data.frame(x2 = sort(runif(n <- 500)))
mymu <- function(x) exp(-2 + 6*sin(2*x-0.2) / (x+0.5)^2)
adata <- transform(adata, y = rpois(n, lambda = mymu(x2)))
mytau <- c(0.25, 0.75); mydof <- 4
fit <- vgam(y ~ s(x2, df = mydof), data=adata, trace=TRUE, maxit = 900,
alaplace2(tau = mytau, llocat = "loglink",
parallel.locat = FALSE))
fitp <- vgam(y ~ s(x2, df = mydof), data = adata, trace=TRUE, maxit=900,
alaplace2(tau = mytau, llocat = "loglink", parallel.locat = TRUE))
par(las = 1); mylwd <- 1.5
with(adata, plot(x2, jitter(y, factor = 0.5), col = "orange",
main = "Example 1; green: parallel.locat = TRUE",
ylab = "y", pch = "o", cex = 0.75))
with(adata, matlines(x2, fitted(fit ), col = "blue",
lty = "solid", lwd = mylwd))
with(adata, matlines(x2, fitted(fitp), col = "green",
lty = "solid", lwd = mylwd))
finexgrid <- seq(0, 1, len = 1001)
for (ii in 1:length(mytau))
lines(finexgrid, qpois(p = mytau[ii], lambda = mymu(finexgrid)),
col = "blue", lwd = mylwd)
fit@extra # Contains useful information
# Example 2: regression quantile at a new tau value from an existing fit
# Nb. regression splines are used here since it is easier.
fitp2 <- vglm(y ~ sm.bs(x2, df = mydof), data = adata, trace = TRUE,
alaplace1(tau = mytau, llocation = "loglink",
parallel.locat = TRUE))
newtau <- 0.5 # Want to refit the model with this tau value
fitp3 <- vglm(y ~ 1 + offset(predict(fitp2)[, 1]),
alaplace1(tau = newtau, llocation = "loglink"), adata)
with(adata, plot(x2, jitter(y, factor = 0.5), col = "orange",
pch = "o", cex = 0.75, ylab = "y",
main = "Example 2; parallel.locat = TRUE"))
with(adata, matlines(x2, fitted(fitp2), col = "blue",
lty = 1, lwd = mylwd))
with(adata, matlines(x2, fitted(fitp3), col = "black",
lty = 1, lwd = mylwd))
# Example 3: noncrossing regression quantiles using a trick: obtain
# successive solutions which are added to previous solutions; use a log
# link to ensure an increasing quantiles at any value of x.
mytau <- seq(0.2, 0.9, by = 0.1)
answer <- matrix(0, nrow(adata), length(mytau)) # Stores the quantiles
adata <- transform(adata, offsety = y*0)
usetau <- mytau
for (ii in 1:length(mytau)) {
# cat("\n\nii = ", ii, "\n")
adata <- transform(adata, usey = y-offsety)
iloc <- ifelse(ii == 1, with(adata, median(y)), 1.0) # Well-chosen!
mydf <- ifelse(ii == 1, 5, 3) # Maybe less smoothing will help
fit3 <- vglm(usey ~ sm.ns(x2, df = mydf), data = adata, trace = TRUE,
alaplace2(tau = usetau[ii], lloc = "loglink", iloc = iloc))
answer[, ii] <- (if(ii == 1) 0 else answer[, ii-1]) + fitted(fit3)
adata <- transform(adata, offsety = answer[, ii])
}
# Plot the results.
with(adata, plot(x2, y, col = "blue",
main = paste("Noncrossing and nonparallel; tau = ",
paste(mytau, collapse = ", "))))
with(adata, matlines(x2, answer, col = "orange", lty = 1))
# Zoom in near the origin.
with(adata, plot(x2, y, col = "blue", xlim = c(0, 0.2), ylim = 0:1,
main = paste("Noncrossing and nonparallel; tau = ",
paste(mytau, collapse = ", "))))
with(adata, matlines(x2, answer, col = "orange", lty = 1))
## End(Not run)
VGAMdata documentation built on Sept. 17, 2024, 9:09 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
View source: R/family.laplace.R
| alaplace | R Documentation |
Asymmetric Laplace Distribution Family Functions
Description
Maximum likelihood estimation of the 1, 2 and 3-parameter asymmetric Laplace distributions (ALDs). The 2-parameter ALD may, with trepidation and lots of skill, sometimes be used as an approximation of quantile regression.
Usage
alaplace1(tau = NULL, llocation = "identitylink",
ilocation = NULL, kappa = sqrt(tau/(1 - tau)), Scale.arg = 1,
ishrinkage = 0.95, parallel.locat = TRUE ~ 0, digt = 4,
idf.mu = 3, zero = NULL, imethod = 1)
alaplace2(tau = NULL, llocation = "identitylink", lscale = "loglink",
ilocation = NULL, iscale = NULL, kappa = sqrt(tau/(1 - tau)),
ishrinkage = 0.95,
parallel.locat = TRUE ~ 0,
parallel.scale = FALSE ~ 0,
digt = 4, idf.mu = 3, imethod = 1, zero = "scale")
alaplace3(llocation = "identitylink", lscale = "loglink",
lkappa = "loglink", ilocation = NULL, iscale = NULL,
ikappa = 1, imethod = 1, zero = c("scale", "kappa"))
Arguments
tau, kappa |
Numeric vectors with
|
llocation, lscale, lkappa |
Character.
Parameter link functions for
location parameter |
ilocation, iscale, ikappa |
Optional initial values. If given, it must be numeric and values are recycled to the appropriate length. The default is to choose the value internally. |
parallel.locat, parallel.scale |
See the |
imethod |
Initialization method. Either the value 1, 2, 3 or 4. |
idf.mu |
Degrees of freedom for the cubic smoothing spline fit
applied to get an initial estimate of the location parameter.
See |
ishrinkage |
How much shrinkage is used when initializing |
Scale.arg |
The value of the scale parameter |
digt |
Passed into |
zero |
See |
Details
These VGAM family functions implement one variant of asymmetric Laplace distributions (ALDs) suitable for quantile regression. Kotz et al. (2001) call it the ALD. Its density function is
f(y;\xi,\sigma,\kappa) = \frac{\sqrt{2}}{\sigma} \,
\frac{\kappa}{1 + \kappa^2} \,
\exp \left( - \frac{\sqrt{2}}{\sigma \, \kappa} |y - \xi |
\right)
for y \leq \xi, and
f(y;\xi,\sigma,\kappa) = \frac{\sqrt{2}}{\sigma} \,
\frac{\kappa}{1 + \kappa^2} \,
\exp \left( - \frac{\sqrt{2} \, \kappa}{\sigma} |y - \xi |
\right)
for y > \xi.
Here, the ranges are for all real y and
\xi, positive \sigma
and positive \kappa.
The special case \kappa = 1 corresponds to the
(symmetric) Laplace distribution of Kotz et al. (2001).
The mean is \xi + \sigma (1/\kappa - \kappa) / \sqrt{2}
and the variance is
\sigma^2 (1 + \kappa^4) / (2 \kappa^2).
The enumeration of the linear/additive predictors used for
alaplace2() is
the first location parameter followed by the first scale parameter,
then the second location parameter followed by the
second scale parameter, etc.
For alaplace3(), only a vector response is handled
and the last (third) linear/additive predictor is for
the asymmetry parameter.
It is known that the maximum likelihood estimate of the
location parameter \xi corresponds to the regression
quantile estimate of the classical quantile regression approach
of Koenker and Bassett (1978). An important property of the
ALD is that
P(Y \leq \xi) = \tau where
\tau = \kappa^2 / (1 + \kappa^2)
so that
\kappa = \sqrt{\tau / (1-\tau)}.
Thus alaplace2() might be used as an alternative to rq
in the quantreg package, although scoring is really
an unsuitable algorithm for estimation here.
Both alaplace1() and alaplace2() can handle
multiple responses, and the number of linear/additive
predictors is dictated by the length of tau or
kappa. The functions alaplace1()
and alaplace2() can also
handle multiple responses (i.e., a matrix response)
but only with a single-valued tau or kappa.
Value
An object of class "vglmff" (see
vglmff-class). The object is used by modelling
functions such as vglm and vgam.
In the extra slot of the fitted object are some list
components which are useful, e.g., the sample proportion of
values which are less than the fitted quantile curves.
Warning
These functions are experimental and especially subject to
change or withdrawal.
The usual MLE regularity conditions do not hold for
this distribution so that misleading inferences may result,
e.g., in the summary and vcov of the object. The
1-parameter ALD can be approximated by extlogF1
which has continuous derivatives and is recommended over
alaplace1.
Care is needed with tau values which are too small,
e.g., for count data with llocation = "loglink" and if
the sample proportion of zeros is greater than tau.
Note
These VGAM family functions use Fisher scoring.
Convergence may be slow and half-stepping is usual
(although one can use trace = TRUE to see which
is the best model and then use maxit to choose
that model)
due to the regularity conditions not holding.
Often the iterations slowly crawl towards the solution so
monitoring the convergence (set trace = TRUE) is highly
recommended.
Instead, extlogF1 is recommended.
For large data sets it is a very good idea to keep the length of
tau/kappa low to avoid large memory requirements.
Then
for parallel.locat = FALSE one can repeatedly fit a model with
alaplace1() with one \tau at a time;
and
for parallel.locat = TRUE one can refit a model with
alaplace1() with one \tau at a time but
using offsets and an intercept-only model.
A second method for solving the noncrossing quantile problem is illustrated below in Example 3. This is called the accumulative quantile method (AQM) and details are in Yee (2015). It does not make the strong parallelism assumption.
The functions alaplace2() and laplace
differ slightly in terms of the parameterizations.
Author(s)
Thomas W. Yee
References
Koenker, R. and Bassett, G. (1978). Regression quantiles. Econometrica, 46, 33–50.
Kotz, S., Kozubowski, T. J. and Podgorski, K. (2001). The Laplace distribution and generalizations: a revisit with applications to communications, economics, engineering, and finance, Boston: Birkhauser.
See Also
ralap,
laplace,
extlogF1,
CommonVGAMffArguments,
lms.bcn,
amlnormal,
sc.studentt2,
simulate.vlm.
Examples
## Not run:
# Example 1: quantile regression with smoothing splines
set.seed(123); adata <- data.frame(x2 = sort(runif(n <- 500)))
mymu <- function(x) exp(-2 + 6*sin(2*x-0.2) / (x+0.5)^2)
adata <- transform(adata, y = rpois(n, lambda = mymu(x2)))
mytau <- c(0.25, 0.75); mydof <- 4
fit <- vgam(y ~ s(x2, df = mydof), data=adata, trace=TRUE, maxit = 900,
alaplace2(tau = mytau, llocat = "loglink",
parallel.locat = FALSE))
fitp <- vgam(y ~ s(x2, df = mydof), data = adata, trace=TRUE, maxit=900,
alaplace2(tau = mytau, llocat = "loglink", parallel.locat = TRUE))
par(las = 1); mylwd <- 1.5
with(adata, plot(x2, jitter(y, factor = 0.5), col = "orange",
main = "Example 1; green: parallel.locat = TRUE",
ylab = "y", pch = "o", cex = 0.75))
with(adata, matlines(x2, fitted(fit ), col = "blue",
lty = "solid", lwd = mylwd))
with(adata, matlines(x2, fitted(fitp), col = "green",
lty = "solid", lwd = mylwd))
finexgrid <- seq(0, 1, len = 1001)
for (ii in 1:length(mytau))
lines(finexgrid, qpois(p = mytau[ii], lambda = mymu(finexgrid)),
col = "blue", lwd = mylwd)
fit@extra # Contains useful information
# Example 2: regression quantile at a new tau value from an existing fit
# Nb. regression splines are used here since it is easier.
fitp2 <- vglm(y ~ sm.bs(x2, df = mydof), data = adata, trace = TRUE,
alaplace1(tau = mytau, llocation = "loglink",
parallel.locat = TRUE))
newtau <- 0.5 # Want to refit the model with this tau value
fitp3 <- vglm(y ~ 1 + offset(predict(fitp2)[, 1]),
alaplace1(tau = newtau, llocation = "loglink"), adata)
with(adata, plot(x2, jitter(y, factor = 0.5), col = "orange",
pch = "o", cex = 0.75, ylab = "y",
main = "Example 2; parallel.locat = TRUE"))
with(adata, matlines(x2, fitted(fitp2), col = "blue",
lty = 1, lwd = mylwd))
with(adata, matlines(x2, fitted(fitp3), col = "black",
lty = 1, lwd = mylwd))
# Example 3: noncrossing regression quantiles using a trick: obtain
# successive solutions which are added to previous solutions; use a log
# link to ensure an increasing quantiles at any value of x.
mytau <- seq(0.2, 0.9, by = 0.1)
answer <- matrix(0, nrow(adata), length(mytau)) # Stores the quantiles
adata <- transform(adata, offsety = y*0)
usetau <- mytau
for (ii in 1:length(mytau)) {
# cat("\n\nii = ", ii, "\n")
adata <- transform(adata, usey = y-offsety)
iloc <- ifelse(ii == 1, with(adata, median(y)), 1.0) # Well-chosen!
mydf <- ifelse(ii == 1, 5, 3) # Maybe less smoothing will help
fit3 <- vglm(usey ~ sm.ns(x2, df = mydf), data = adata, trace = TRUE,
alaplace2(tau = usetau[ii], lloc = "loglink", iloc = iloc))
answer[, ii] <- (if(ii == 1) 0 else answer[, ii-1]) + fitted(fit3)
adata <- transform(adata, offsety = answer[, ii])
}
# Plot the results.
with(adata, plot(x2, y, col = "blue",
main = paste("Noncrossing and nonparallel; tau = ",
paste(mytau, collapse = ", "))))
with(adata, matlines(x2, answer, col = "orange", lty = 1))
# Zoom in near the origin.
with(adata, plot(x2, y, col = "blue", xlim = c(0, 0.2), ylim = 0:1,
main = paste("Noncrossing and nonparallel; tau = ",
paste(mytau, collapse = ", "))))
with(adata, matlines(x2, answer, col = "orange", lty = 1))
## End(Not run)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.