fitGEV: Fit a Generalized Extreme value (GEV) GAMLSS Model

Description

Fits a Generalized Additive Model (GAM) for Location, Scale and Shape with a GEV response distribution, using the function gamlss::gamlss().

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,
  ...
)

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").

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.

Examples

Run this code

# Load gamlss, for the functions term.plot() and wp()
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 ~ 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)