loglaplace: Log-Laplace and Logit-Laplace Distribution Family Functions

Description Usage Arguments Details Value Warning Note Author(s) References See Also Examples

Description

Maximum likelihood estimation of the 1-parameter log-Laplace and the 1-parameter logit-Laplace distributions. These may be used for quantile regression for counts and proportions respectively.

Usage

1
2
3
4
5
6
7
8
9
loglaplace1(tau = NULL, llocation = "loglink",
    ilocation = NULL, kappa = sqrt(tau/(1 - tau)), Scale.arg = 1,
    ishrinkage = 0.95, parallel.locat = FALSE, digt = 4,
    idf.mu = 3, rep0 = 0.5, minquantile = 0, maxquantile = Inf,
    imethod = 1, zero = NULL)
logitlaplace1(tau = NULL, llocation = "logitlink",
    ilocation = NULL, kappa = sqrt(tau/(1 - tau)),
    Scale.arg = 1, ishrinkage = 0.95, parallel.locat = FALSE,
    digt = 4, idf.mu = 3, rep01 = 0.5, imethod = 1, zero = NULL)

Arguments

tau, kappa

See alaplace1.

llocation

Character. Parameter link functions for location parameter xi. See Links for more choices. However, this argument should be left unchanged with count data because it restricts the quantiles to be positive. With proportions data llocation can be assigned a link such as logitlink, probitlink, clogloglink, etc.

ilocation

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

Logical. Should the quantiles be parallel on the transformed scale (argument llocation)? Assigning this argument to TRUE circumvents the seriously embarrassing quantile crossing problem.

imethod

Initialization method. Either the value 1, 2, or ....

idf.mu, ishrinkage, Scale.arg, digt, zero

See alaplace1.

rep0, rep01

Numeric, positive. Replacement values for 0s and 1s respectively. For count data, values of the response whose value is 0 are replaced by rep0; it avoids computing log(0). For proportions data values of the response whose value is 0 or 1 are replaced by min(rangey01[1]/2, rep01/w[y< = 0]) and max((1 + rangey01[2])/2, 1-rep01/w[y >= 1]) respectively; e.g., it avoids computing logitlink(0) or logitlink(1). Here, rangey01 is the 2-vector range(y[(y > 0) & (y < 1)]) of the response.

minquantile, maxquantile

Numeric. The minimum and maximum values possible in the quantiles. These argument are effectively ignored by default since loglink keeps all quantiles positive. However, if llocation = logofflink(offset = 1) then it is possible that the fitted quantiles have value 0 because minquantile = 0.

Details

These VGAM family functions implement translations of the asymmetric Laplace distribution (ALD). The resulting variants may be suitable for quantile regression for count data or sample proportions. For example, a log link applied to count data is assumed to follow an ALD. Another example is a logit link applied to proportions data so as to follow an ALD. A positive random variable Y is said to have a log-Laplace distribution if Y = exp(W) where W has an ALD. There are many variants of ALDs and the one used here is described in alaplace1.

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. For example, the sample proportion of values which are less than the fitted quantile curves, which is sum(wprior[y <= location]) / sum(wprior) internally. Here, wprior are the prior weights (called ssize below), y is the response and location is a fitted quantile curve. This definition comes about naturally from the transformed ALD data.

Warning

The VGAM family function logitlaplace1 will not handle a vector of just 0s and 1s as the response; it will only work satisfactorily if the number of trials is large.

See alaplace1 for other warnings. Care is needed with tau values which are too small, e.g., for count data the sample proportion of zeros must be less than all values in tau. Similarly, this also holds with logitlaplace1, which also requires all tau values to be less than the sample proportion of ones.

Note

The form of input for logitlaplace1 as response is a vector of proportions (values in [0,1]) and the number of trials is entered into the weights argument of vglm/vgam. See Example 2 below. See alaplace1 for other notes in general.

Author(s)

Thomas W. Yee

References

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.

Kozubowski, T. J. and Podgorski, K. (2003). Log-Laplace distributions. International Mathematical Journal, 3, 467–495.

Yee, T. W. (2020). Quantile regression for counts and proportions. In preparation.

See Also

alaplace1, dloglap.

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
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
# Example 1: quantile regression of counts with regression splines
set.seed(123); my.k <- exp(0)
adata <- data.frame(x2 = sort(runif(n <- 500)))
mymu <- function(x) exp( 1 + 3*sin(2*x) / (x+0.5)^2)
adata <- transform(adata, y = rnbinom(n, mu = mymu(x2), size = my.k))
mytau <- c(0.1, 0.25, 0.5, 0.75, 0.9); mydof = 3
# halfstepping is usual:
fitp <- vglm(y ~ sm.bs(x2, df = mydof), data = adata, trace = TRUE,
            loglaplace1(tau = mytau, parallel.locat = TRUE))

## Not run: 
par(las = 1)  # Plot on a log1p() scale
mylwd <- 1.5


plot(jitter(log1p(y), factor = 1.5) ~ x2, adata, col = "red", pch = "o",
     main = "Example 1; darkgreen=truth, blue=estimated", cex = 0.75)
with(adata, matlines(x2, log1p(fitted(fitp)), col = "blue",
                      lty = 1, lwd = mylwd))
finexgrid <- seq(0, 1, len = 201)
for (ii in 1:length(mytau))
  lines(finexgrid, col = "darkgreen", lwd = mylwd,
        log1p(qnbinom(p = mytau[ii], mu = mymu(finexgrid), si = my.k)))

## End(Not run)
fitp@extra  # Contains useful information


# Example 2: sample proportions
set.seed(123); nnn <- 1000; ssize <- 100  # ssize = 1 will not work!
adata <- data.frame(x2 = sort(runif(nnn)))
mymu <- function(x) logitlink( 1.0 + 4*x, inv = TRUE)
adata <- transform(adata, ssize = ssize,
                    y2 = rbinom(nnn, size = ssize, prob = mymu(x2)) / ssize)

mytau <- c(0.25, 0.50, 0.75)
fit1 <- vglm(y2 ~ sm.bs(x2, df = 3), data = adata,
             weights = ssize, trace = TRUE,
             logitlaplace1(tau = mytau, lloc = "clogloglink", paral = TRUE))

## Not run: 
# Check the solution.  Note: this may be like comparing apples with oranges.
plotvgam(fit1, se = TRUE, scol = "red", lcol = "blue",
         main = "Truth = 'darkgreen'")
# Centered approximately !
linkFunctionChar <- as.character(fit1@misc$link)
adata <- transform(adata, trueFunction=
                   theta2eta(theta = mymu(x2), link=linkFunctionChar))
with(adata, lines(x2, trueFunction - mean(trueFunction), col = "darkgreen"))


# Plot the data + fitted quantiles (on the original scale)
myylim <- with(adata, range(y2))
plot(y2 ~ x2, adata, col = "blue", ylim = myylim, las = 1,
     pch = ".", cex = 2.5)
with(adata, matplot(x2, fitted(fit1), add = TRUE, lwd = 3, type = "l"))
truecol <- rep(1:3, len = fit1@misc$M)  # Add the 'truth'
smallxgrid <- seq(0, 1, len = 501)
for (ii in 1:length(mytau))
  lines(smallxgrid, col = truecol[ii], lwd = 2,
        qbinom(p = mytau[ii], prob = mymu(smallxgrid), size = ssize) / ssize)


# Plot on the eta (== logitlink()/probit()/...) scale
with(adata, matplot(x2, predict(fit1), add = FALSE, lwd = 3, type = "l"))
# Add the 'truth'
for (ii in 1:length(mytau)) {
  true.quant <- qbinom(p = mytau[ii], pr = mymu(smallxgrid), si = ssize) / ssize
  lines(smallxgrid, theta2eta(theta = true.quant, link = linkFunctionChar),
        col = truecol[ii], lwd = 2)
}

## End(Not run)

Example output

Loading required package: stats4
Loading required package: splines
VGLM    linear loop  1 :  loglikelihood = -4627.9133
VGLM    linear loop  2 :  loglikelihood = -4567.4494
VGLM    linear loop  3 :  loglikelihood = -4543.2723
VGLM    linear loop  4 :  loglikelihood = -4533.293
VGLM    linear loop  5 :  loglikelihood = -4529.5149
VGLM    linear loop  6 :  loglikelihood = -4527.8516
VGLM    linear loop  7 :  loglikelihood = -4527.3429
VGLM    linear loop  8 :  loglikelihood = -4527.3436
Taking a modified step.
VGLM    linear loop  8 :  loglikelihood = -4527.2428
VGLM    linear loop  9 :  loglikelihood = -4527.1886
VGLM    linear loop  10 :  loglikelihood = -4527.1903
Taking a modified step.
VGLM    linear loop  10 :  loglikelihood = -4526.952
VGLM    linear loop  11 :  loglikelihood = -4527.1339
Taking a modified step..
VGLM    linear loop  11 :  loglikelihood = -4526.9166
VGLM    linear loop  12 :  loglikelihood = -4527.2443
Taking a modified step.....
VGLM    linear loop  12 :  loglikelihood = -4526.9138
VGLM    linear loop  13 :  loglikelihood = -4527.0252
Taking a modified step....
VGLM    linear loop  13 :  loglikelihood = -4526.9119
VGLM    linear loop  14 :  loglikelihood = -4526.9658
Taking a modified step......
VGLM    linear loop  14 :  loglikelihood = -4526.9118
VGLM    linear loop  15 :  loglikelihood = -4527.026
Taking a modified step......
VGLM    linear loop  15 :  loglikelihood = -4526.9118
Warning message:
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
$M
[1] 5

$Scale
[1] 1 1 1 1 1

$kappa
[1] 0.3333333 0.5773503 1.0000000 1.7320508 3.0000000

$tau
[1] 0.10 0.25 0.50 0.75 0.90

$n
[1] 500

$y.names
[1] "tau = 0.1"  "tau = 0.25" "tau = 0.5"  "tau = 0.75" "tau = 0.9" 

$individual
[1] FALSE

$Scale.arg
[1] 1

$percentile
[1] 10.6 25.0 50.0 75.0 90.0

VGLM    linear loop  1 :  loglikelihood = -209353.75
VGLM    linear loop  2 :  loglikelihood = -197285.15
VGLM    linear loop  3 :  loglikelihood = -241183.1
Taking a modified step.
VGLM    linear loop  3 :  loglikelihood = -191845.6
VGLM    linear loop  4 :  loglikelihood = -213969.07
Taking a modified step.
VGLM    linear loop  4 :  loglikelihood = -181088.15
VGLM    linear loop  5 :  loglikelihood = -192869.9
Taking a modified step..
VGLM    linear loop  5 :  loglikelihood = -179741.64
VGLM    linear loop  6 :  loglikelihood = -180247.84
Taking a modified step.
VGLM    linear loop  6 :  loglikelihood = -179681.87
VGLM    linear loop  7 :  loglikelihood = -179945.99
Taking a modified step..
VGLM    linear loop  7 :  loglikelihood = -179658.56
VGLM    linear loop  8 :  loglikelihood = -179666.97
Taking a modified step.
VGLM    linear loop  8 :  loglikelihood = -179656.91
VGLM    linear loop  9 :  loglikelihood = -179662.7
Taking a modified step..
VGLM    linear loop  9 :  loglikelihood = -179656.74
VGLM    linear loop  10 :  loglikelihood = -179657.69
Taking a modified step..
VGLM    linear loop  10 :  loglikelihood = -179656.52
VGLM    linear loop  11 :  loglikelihood = -179658.8
Taking a modified step...
VGLM    linear loop  11 :  loglikelihood = -179656.41
VGLM    linear loop  12 :  loglikelihood = -179658.66
Taking a modified step....
VGLM    linear loop  12 :  loglikelihood = -179656.39
VGLM    linear loop  13 :  loglikelihood = -179658.59
Taking a modified step....
VGLM    linear loop  13 :  loglikelihood = -179656.38
VGLM    linear loop  14 :  loglikelihood = -179657.08
Taking a modified step...
VGLM    linear loop  14 :  loglikelihood = -179656.37
VGLM    linear loop  15 :  loglikelihood = -179659.31
Taking a modified step....
VGLM    linear loop  15 :  loglikelihood = -179656.35
VGLM    linear loop  16 :  loglikelihood = -179658.7
Taking a modified step......
VGLM    linear loop  16 :  loglikelihood = -179656.34
VGLM    linear loop  17 :  loglikelihood = -179658.73
Taking a modified step.....
VGLM    linear loop  17 :  loglikelihood = -179656.34
VGLM    linear loop  18 :  loglikelihood = -179658.62
Taking a modified step......
VGLM    linear loop  18 :  loglikelihood = -179656.34
VGLM    linear loop  19 :  loglikelihood = -179656.85
Taking a modified step....
VGLM    linear loop  19 :  loglikelihood = -179656.33
VGLM    linear loop  20 :  loglikelihood = -179658.98
Taking a modified step......
VGLM    linear loop  20 :  loglikelihood = -179656.32
VGLM    linear loop  21 :  loglikelihood = -179656.86
Taking a modified step.....
VGLM    linear loop  21 :  loglikelihood = -179656.32
VGLM    linear loop  22 :  loglikelihood = -179657.13
Taking a modified step......
VGLM    linear loop  22 :  loglikelihood = -179656.32
VGLM    linear loop  23 :  loglikelihood = -179658.63
Taking a modified step........
VGLM    linear loop  23 :  loglikelihood = -179656.31
VGLM    linear loop  24 :  loglikelihood = -179657.53
Taking a modified step........
VGLM    linear loop  24 :  loglikelihood = -179656.31
VGLM    linear loop  25 :  loglikelihood = -179656.87
Taking a modified step...............
VGLM    linear loop  25 :  loglikelihood = -179656.31
Warning message:
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 documentation built on Jan. 16, 2021, 5:21 p.m.