fitGEV {gamlssx}

Fit a Generalized Extreme value (GEV) GAMLSS Model

Description

Describe

Usage

fitGEV(
  formula,
  data,
  scoring = c("fisher", "quasi"),
  mu.link = "identity",
  sigma.link = "log",
  xi.link = "identity",
  stepLength = 1,
  stepAttempts = 2,
  stepReduce = 2,
  steps = FALSE,
  ...
)

Arguments

formula

a formula object, with the response on the left of an ~ operator, and the terms, separated by + operators, on the right. Nonparametric smoothing terms are indicated by pb() for penalised beta splines, cs for smoothing splines, lo for loess smooth terms and random or ra for random terms, e.g. y~cs(x,df=5)+x1+x2*x3. Additional smoothers can be added by creating the appropriate interface. Interactions with nonparametric smooth terms are not fully supported, but will not produce errors; they will simply produce the usual parametric interaction

data

a data frame containing the variables occurring in the formula, e.g. data=aids. If this is missing, the variables should be on the search list.

scoring

A character scalar. If scoring = "fisher" then the weights used in the fitting algorithm are based on the expected Fisher information, that is, a Fisher's scoring algorithm is used. If scoring = "quasi" then these weights are based on the cross products of the first derivatives of the log-likelihood, leading to a quasi Newton scoring algorithm.

mu.link, sigma.link, xi.link

Character scalars to set the respective link functions for the location (mu), scale (sigma) and shape (xi) parameters. The latter is passed to gamlss::gamlss() as nu.link.

stepLength

A numeric vector of positive values. The initial values of the step lengths mu.step, sigma.step and nu.step passed to gamlss::gamlss.control() in the first attempt to fit the model by calling gamlss::gamlss(). If stepLength has a length that is less than 3 then stepLength is recycled to have length 3.

stepAttempts

A non-negative integer. If the first call to gamlss::gamlss() throws an error then we make stepAttempts further attempts to fit the model, each time dividing by 2 the values of mu.step, sigma.step and nu.step supplied to gamlss::gamlss.control(). If stepAttempts < 1 then no further attempts are made.

stepReduce

A number greater than 1. The factor by which the step lengths in stepLength are reduced for each extra attempt to fit the model. The default, stepReduce = 2 means that the step lengths are halved for each extra attempt.

steps

A logical scalar. Pass steps = TRUE to write to the console the current value of stepLength for each call to gamlss::gamlss().

...

Further arguments passed to gamlss::gamlss(), in particular method, which sets the fitting algorithm, with options RS(), CG() or mixed(). The default, method = RS() seems to work well, as does method = mixed(). In contrast, method = CG() often requires the step length to be reduced before convergence is achieved. fitGEV() attempts to do this automatically. See stepAttempts. Pass trace = FALSE (to gamlss::gamlss.control()) to avoid writing to the console the global deviance after each outer iteration of the gamlss fitting algorithm.

Details

See gamlss::gamlss() for information about the model and the fitting algorithm.

Value

Returns a gamlss object. See the Value section of gamlss::gamlss(). The class of the returned object is c("gamlssx", "gamlss", "gam", "glm", "lm").

See Also

GEV, gamlss.dist::gamlss.family(), gamlss::gamlss()

Examples

# Load gamlss, for the function pb()
library(gamlss)

##### Simulated data

set.seed(17012023)
n <- 100
x <- stats::runif(n)
mu <- 1 + 2 * x
sigma <- 1
xi <- 0.25
y <- nieve::rGEV(n = 1, loc = mu, scale = sigma, shape = xi)
data <- data.frame(y = as.numeric(y), x = x)
plot(x, y)

# Fit model using the default RS method with Fisher's scoring
mod <- fitGEV(y ~ gamlss::pb(x), data = data)
# Summary of model fit
summary(mod)
# Residual diagnostic plots
plot(mod, xlab = "x", ylab = "y")
# Data plus fitted curve
plot(data$x, data$y, xlab = "x", ylab = "y")
lines(data$x, fitted(mod))

# Fit model using the mixed method and quasi-Newton scoring
# Use trace = FALSE to prevent writing the global deviance to the console
mod <- fitGEV(y ~ pb(x), data = data, method = mixed(), scoring = "quasi",
              trace = FALSE)

# Fit model using the CG method
# The default step length of 1 needs to be reduced to enable convergence
# Use steps = TRUE to write the step lengths to the console
mod <- fitGEV(y ~ pb(x), data = data, method = CG(), steps = TRUE)

##### Fremantle annual maximum sea levels
##### See also the gamlssx package README file

# Transform Year so that it is centred on 0
fremantle <- transform(fremantle, cYear = Year - median(Year))

# Plot sea level against year and against SOI
plot(fremantle$Year, fremantle$SeaLevel, xlab = "year", ylab = "sea level (m)")
plot(fremantle$SOI, fremantle$SeaLevel, xlab = "SOI", ylab = "sea level (m)")

# Fit a model with P-spline effects of cYear and SOI on location and scale
# The default links are identity for location and log for scale
mod <- fitGEV(SeaLevel ~ pb(cYear) + pb(SOI),
             sigma.formula = ~ pb(cYear) + pb(SOI),
             data = fremantle)

# Summary of model fit
summary(mod)
# Model diagnostic plots
plot(mod)
# Worm plot
wp(mod)
# Plot of the fitted component smooth functions
# Note: gamlss::term.plot() does not include uncertainty about the intercept
# Location mu
term.plot(mod, rug = TRUE, pages = 1)
# Scale sigma
term.plot(mod, what = "sigma", rug = TRUE, pages = 1)