alaplace: 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 $\tau$th 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.

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.

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.

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.

Examples

Run this code

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

Run the code above in your browser using DataLab