fgev: Maximum-likelihood Fitting of the Generalized Extreme Value...

Description Usage Arguments Details Value Warning References See Also Examples

Description

Maximum-likelihood fitting for the generalized extreme value distribution, including linear modelling of the location parameter, and allowing any of the parameters to be held fixed if desired.

Usage

1
2
fgev(x, start, ..., nsloc = NULL, prob = NULL, std.err = TRUE,
    corr = FALSE, method = "BFGS", warn.inf = TRUE)

Arguments

x

A numeric vector, which may contain missing values.

start

A named list giving the initial values for the parameters over which the likelihood is to be maximized. If start is omitted the routine attempts to find good starting values using moment estimators.

...

Additional parameters, either for the GEV model or for the optimization function optim. If parameters of the model are included they will be held fixed at the values given (see Examples).

nsloc

A data frame with the same number of rows as the length of x, for linear modelling of the location parameter. The data frame is treated as a covariate matrix (excluding the intercept). A numeric vector can be given as an alternative to a single column data frame.

prob

Controls the parameterization of the model (see Details). Should be either NULL (the default), or a probability in the closed interval [0,1].

std.err

Logical; if TRUE (the default), the standard errors are returned.

corr

Logical; if TRUE, the correlation matrix is returned.

method

The optimization method (see optim for details).

warn.inf

Logical; if TRUE (the default), a warning is given if the negative log-likelihood is infinite when evaluated at the starting values.

Details

If prob is NULL (the default):

For stationary models the parameter names are loc, scale and shape, for the location, scale and shape parameters respectively. For non-stationary models, the parameter names are loc, locx1, ..., locxn, scale and shape, where x1, ..., xn are the column names of nsloc, so that loc is the intercept of the linear model, and locx1, ..., locxn are the ncol(nsloc) coefficients. If nsloc is a vector it is converted into a single column data frame with column name trend, and hence the associated trend parameter is named loctrend.

If \code{prob} = p is a probability:

The fit is performed using a different parameterization. Let a, b and s denote the location, scale and shape parameters of the GEV distribution. For stationary models, the distribution is parameterized using (z_p,b,s), where

z_p = a - b/s (1 - (-log(1 - p))^s)

is such that G(z_p) = 1 - p, where G is the GEV distribution function. \code{prob} = p is therefore the probability in the upper tail corresponding to the quantile z_p. If prob is zero, then z_p is the upper end point a - b/s, and s is restricted to the negative (Weibull) axis. If prob is one, then z_p is the lower end point a - b/s, and s is restricted to the positive (Frechet) axis. The parameter names are quantile, scale and shape, for z_p, b and s respectively.

For non-stationary models the parameter z_p is again given by the equation above, but a becomes the intercept of the linear model for the location parameter, so that quantile replaces (the intercept) loc, and hence the parameter names are quantile, locx1, ..., locxn, scale and shape, where x1, ..., xn are the column names of nsloc.

In either case:

For non-stationary fitting it is recommended that the covariates within the linear model for the location parameter are (at least approximately) centered and scaled (i.e.\ that the columns of nsloc are centered and scaled), particularly if automatic starting values are used, since the starting values for the associated parameters are then zero.

Value

Returns an object of class c("gev","uvevd","evd").

The generic accessor functions fitted (or fitted.values), std.errors, deviance, logLik and AIC extract various features of the returned object.

The functions profile and profile2d are used to obtain deviance profiles for the model parameters. In particular, profiles of the quantile z_p can be calculated and plotted when \code{prob} = p. The function anova compares nested models. The function plot produces diagnostic plots.

An object of class c("gev","uvevd","evd") is a list containing at most the following components

estimate

A vector containing the maximum likelihood estimates.

std.err

A vector containing the standard errors.

fixed

A vector containing the parameters of the model that have been held fixed.

param

A vector containing all parameters (optimized and fixed).

deviance

The deviance at the maximum likelihood estimates.

corr

The correlation matrix.

var.cov

The variance covariance matrix.

convergence, counts, message

Components taken from the list returned by optim.

data

The data passed to the argument x.

tdata

The data, transformed to stationarity (for non-stationary models).

nsloc

The argument nsloc.

n

The length of x.

prob

The argument prob.

loc

The location parameter. If prob is NULL (the default), this will also be an element of param.

call

The call of the current function.

Warning

The standard errors and the correlation matrix in the returned object are taken from the observed information, calculated by a numerical approximation. They must be interpreted with caution when the shape parameter is less than -0.5, because the usual asymptotic properties of maximum likelihood estimators do not then hold (Smith, 1985).

References

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

See Also

anova.evd, optim, plot.uvevd, profile.evd, profile2d.evd

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
uvdata <- rgev(100, loc = 0.13, scale = 1.1, shape = 0.2)
trend <- (-49:50)/100
M1 <- fgev(uvdata, nsloc = trend, control = list(trace = 1))
M2 <- fgev(uvdata)
M3 <- fgev(uvdata, shape = 0)
M4 <- fgev(uvdata, scale = 1, shape = 0)
anova(M1, M2, M3, M4)
par(mfrow = c(2,2))
plot(M2)
## Not run: M2P <- profile(M2)
## Not run: plot(M2P)

rnd <- runif(100, min = -.5, max = .5)
fgev(uvdata, nsloc = data.frame(trend = trend, random = rnd))
fgev(uvdata, nsloc = data.frame(trend = trend, random = rnd), locrandom = 0)

uvdata <- rgev(100, loc = 0.13, scale = 1.1, shape = 0.2)
M1 <- fgev(uvdata, prob = 0.1)
M2 <- fgev(uvdata, prob = 0.01)
## Not run: M1P <- profile(M1, which = "quantile")
## Not run: M2P <- profile(M2, which = "quantile")
## Not run: plot(M1P)
## Not run: plot(M2P)

Example output

initial  value 212.891896 
iter  10 value 191.536175
iter  10 value 191.536173
final  value 191.536173 
converged
Analysis of Deviance Table

   M.Df Deviance Df   Chisq Pr(>chisq)    
M1    4   383.07                          
M2    3   384.39  1  1.3159     0.2513    
M3    2   401.38  1 16.9914  3.755e-05 ***
M4    1   423.26  1 21.8783  2.905e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Call: fgev(x = uvdata, nsloc = data.frame(trend = trend, random = rnd)) 
Deviance: 381.8589 

Estimates
      loc   loctrend  locrandom      scale      shape  
   0.2507    -0.5075     0.3816     1.1768     0.2985  

Standard Errors
      loc   loctrend  locrandom      scale      shape  
  0.14031    0.37801    0.33744    0.11663    0.09975  

Optimization Information
  Convergence: successful 
  Function Evaluations: 28 
  Gradient Evaluations: 12 


Call: fgev(x = uvdata, locrandom = 0, nsloc = data.frame(trend = trend,      random = rnd)) 
Deviance: 383.0723 

Estimates
     loc  loctrend     scale     shape  
  0.2852   -0.4501    1.1963    0.2794  

Standard Errors
     loc  loctrend     scale     shape  
 0.13851   0.39379   0.11568   0.09341  

Optimization Information
  Convergence: successful 
  Function Evaluations: 25 
  Gradient Evaluations: 10 

evd documentation built on April 25, 2018, 5:04 p.m.