Learn R Programming

VGAM (version 1.1-2)

gev: Generalized Extreme Value Regression Family Function

Description

Maximum likelihood estimation of the 3-parameter generalized extreme value (GEV) distribution.

Usage

gev(llocation = "identitylink", lscale = "loglink",
    lshape = logofflink(offset = 0.5), percentiles = c(95, 99),
    ilocation = NULL, iscale = NULL, ishape = NULL, imethod = 1,
    gprobs.y = (1:9)/10, gscale.mux = exp((-5:5)/6),
    gshape = (-5:5) / 11 + 0.01,
    iprobs.y = NULL, tolshape0 = 0.001,
    type.fitted = c("percentiles", "mean"),
    zero = c("scale", "shape"))
gevff(llocation = "identitylink", lscale = "loglink",
    lshape = logofflink(offset = 0.5), percentiles = c(95, 99),
    ilocation = NULL, iscale = NULL, ishape = NULL, imethod = 1,
    gprobs.y = (1:9)/10, gscale.mux = exp((-5:5)/6),
    gshape = (-5:5) / 11 + 0.01,
    iprobs.y = NULL, tolshape0 = 0.001,
    type.fitted = c("percentiles", "mean"), zero = c("scale", "shape"))

Arguments

llocation, lscale, lshape

Parameter link functions for \(\mu\), \(\sigma\) and \(\xi\) respectively. See Links for more choices.

For the shape parameter, the default logofflink link has an offset called \(A\) below; and then the linear/additive predictor is \(\log(\xi+A)\) which means that \(\xi > -A\). For technical reasons (see Details) it is a good idea for \(A = 0.5\).

percentiles

Numeric vector of percentiles used for the fitted values. Values should be between 0 and 100. This argument is ignored if type.fitted = "mean".

type.fitted

See CommonVGAMffArguments for information. The default is to use the percentiles argument. If "mean" is chosen, then the mean \(\mu + \sigma (\Gamma(1-\xi)-1) / \xi\) is returned as the fitted values, and these are only defined for \(\xi<1\).

ilocation, iscale, ishape

Numeric. Initial value for the location parameter, \(\sigma\) and \(\xi\). A NULL means a value is computed internally. The argument ishape is more important than the other two. If a failure to converge occurs, or even to obtain initial values occurs, try assigning ishape some value (positive or negative; the sign can be very important). Also, in general, a larger value of iscale tends to be better than a smaller value.

imethod

Initialization method. Either the value 1 or 2. If both methods fail then try using ishape. See CommonVGAMffArguments for information.

gshape

Numeric vector. The values are used for a grid search for an initial value for \(\xi\). See CommonVGAMffArguments for information.

gprobs.y, gscale.mux, iprobs.y

Numeric vectors, used for the initial values. See CommonVGAMffArguments for information.

tolshape0

Passed into dgev when computing the log-likelihood.

zero

A specifying which linear/additive predictors are modelled as intercepts only. The values can be from the set {1,2,3} corresponding respectively to \(\mu\), \(\sigma\), \(\xi\). If zero = NULL then all linear/additive predictors are modelled as a linear combination of the explanatory variables. For many data sets having zero = 3 is a good idea. See CommonVGAMffArguments for information.

Value

An object of class "vglmff" (see vglmff-class). The object is used by modelling functions such as vglm, and vgam.

Warning

Currently, if an estimate of \(\xi\) is too close to 0 then an error may occur for gev() with multivariate responses. In general, gevff() is more reliable than gev().

Fitting the GEV by maximum likelihood estimation can be numerically fraught. If \(1 + \xi (y-\mu)/ \sigma \leq 0\) then some crude evasive action is taken but the estimation process can still fail. This is particularly the case if vgam with s is used; then smoothing is best done with vglm with regression splines (bs or ns) because vglm implements half-stepsizing whereas vgam doesn't (half-stepsizing helps handle the problem of straying outside the parameter space.)

Details

The GEV distribution function can be written $$G(y) = \exp( -[ (y-\mu)/ \sigma ]_{+}^{- 1/ \xi}) $$ where \(\sigma > 0\), \(-\infty < \mu < \infty\), and \(1 + \xi(y-\mu)/\sigma > 0\). Here, \(x_+ = \max(x,0)\). The \(\mu\), \(\sigma\), \(\xi\) are known as the location, scale and shape parameters respectively. The cases \(\xi>0\), \(\xi<0\), \(\xi = 0\) correspond to the Frechet, reverse Weibull, and Gumbel types respectively. It can be noted that the Gumbel (or Type I) distribution accommodates many commonly-used distributions such as the normal, lognormal, logistic, gamma, exponential and Weibull.

For the GEV distribution, the \(k\)th moment about the mean exists if \(\xi < 1/k\). Provided they exist, the mean and variance are given by \(\mu+\sigma\{ \Gamma(1-\xi)-1\}/ \xi\) and \(\sigma^2 \{ \Gamma(1-2\xi) - \Gamma^2(1-\xi) \} / \xi^2\) respectively, where \(\Gamma\) is the gamma function.

Smith (1985) established that when \(\xi > -0.5\), the maximum likelihood estimators are completely regular. To have some control over the estimated \(\xi\) try using lshape = logofflink(offset = 0.5), say, or lshape = extlogitlink(min = -0.5, max = 0.5), say.

References

Yee, T. W. and Stephenson, A. G. (2007) Vector generalized linear and additive extreme value models. Extremes, 10, 1--19.

Tawn, J. A. (1988) An extreme-value theory model for dependent observations. Journal of Hydrology, 101, 227--250.

Prescott, P. and Walden, A. T. (1980) Maximum likelihood estimation of the parameters of the generalized extreme-value distribution. Biometrika, 67, 723--724.

Smith, R. L. (1985) Maximum likelihood estimation in a class of nonregular cases. Biometrika, 72, 67--90.

See Also

rgev, gumbel, gumbelff, guplot, rlplot.gevff, gpd, weibullR, frechet, extlogitlink, oxtemp, venice, CommonVGAMffArguments.

Examples

Run this code
# NOT RUN {
# Multivariate example
fit1 <- vgam(cbind(r1, r2) ~ s(year, df = 3), gev(zero = 2:3),
             data = venice, trace = TRUE)
coef(fit1, matrix = TRUE)
head(fitted(fit1))
par(mfrow = c(1, 2), las = 1)
plot(fit1, se = TRUE, lcol = "blue", scol = "forestgreen",
     main = "Fitted mu(year) function (centered)", cex.main = 0.8)
with(venice, matplot(year, depvar(fit1)[, 1:2], ylab = "Sea level (cm)",
     col = 1:2, main = "Highest 2 annual sea levels", cex.main = 0.8))
with(venice, lines(year, fitted(fit1)[,1], lty = "dashed", col = "blue"))
legend("topleft", lty = "dashed", col = "blue", "Fitted 95 percentile")

# Univariate example
(fit <- vglm(maxtemp ~ 1, gevff, data = oxtemp, trace = TRUE))
head(fitted(fit))
coef(fit, matrix = TRUE)
Coef(fit)
vcov(fit)
vcov(fit, untransform = TRUE)
sqrt(diag(vcov(fit)))  # Approximate standard errors
rlplot(fit)
# }

Run the code above in your browser using DataLab