gamlssx

Generalised Additive Extreme Value Models for Location, Scale and Shape

https://github.com/paulnorthrop/gamlssx

Keywords

Repository

Generalised Additive Extreme Value Models for Location, Scale and Shape

README.Rmd

---
output: github_document
---



```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "man/figures/README-",
  out.width = "100%"
)
```

# gamlssx 

[![Appveyor Build status](https://ci.appveyor.com/api/projects/status/99jojhgk9t4agdmv/branch/main?svg=true)](https://ci.appveyor.com/project/paulnorthrop/gamlssx/branch/main)
[![R-CMD-check](https://github.com/paulnorthrop/gamlssx/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/paulnorthrop/gamlssx/actions/workflows/R-CMD-check.yaml)
[![Coverage Status](https://codecov.io/github/paulnorthrop/gamlssx/coverage.svg?branch=master)](https://app.codecov.io/github/paulnorthrop/gamlssx?branch=master)
[![CRAN_Status_Badge](https://www.r-pkg.org/badges/version/gamlssx)](https://cran.r-project.org/package=gamlssx)
[![Downloads (monthly)](https://cranlogs.r-pkg.org/badges/gamlssx?color=brightgreen)](https://cran.r-project.org/package=gamlssx)
[![Downloads (total)](https://cranlogs.r-pkg.org/badges/grand-total/gamlssx?color=brightgreen)](https://cran.r-project.org/package=gamlssx)

## Generalized Additive Extreme Value Models for Location, Scale and Shape

The main aim of the `gamlssx` package is to enable a generalized extreme value (GEV) to be used as the response distribution in a generalized additive model for location scale and shape (GAMLSS), as implemented in the [gamlss](https://CRAN.R-project.org/package=gamlss) R package. The [gamlss.dist](https://CRAN.R-project.org/package=gamlss.dist) R package does offer reversed GEV distribution via in `RGE` family, but (a) this is not the usual parameterization of a GEV distribution (for block maxima), and (b) in `RGE`, the shape parameter is restricted to have a particular sign, which is undesirable because the sign of the shape parameter influences strongly extremal behaviour. The `gamlssx` package uses the usual parameterization, with a shape parameter $\xi$, and imposes only the restriction that, for each observation in the data, $\xi > -1/2$, which is necessary for the usual asymptotic likelihood theory to be applicable.

See [Rigby and Stasinopoulos (2005)](https://doi.org/10.1111%2Fj.1467-9876.2005.00510.x) and the [gamlss home page](https://www.gamlss.com/) for details of the GAMLSS methodology. See also Gavin Simpson's blog post [Modelling extremes using generalized additive models](https://fromthebottomoftheheap.net/2017/01/25/modelling-extremes-with-gams/) for an overview of the use of GAMs for modelling extreme values, which uses the [mgcv](https://cran.r-project.org/package=mgcv) R package to fit similar models. The [VGAM](https://CRAN.R-project.org/package=VGAM) and [evgam](https://CRAN.R-project.org/package=evgam) R packages can also be used.

## An example

We consider the `fremantle` data include in the `gamlssx` package, which is a copy of data of the same name from the [`ismev` R package](https://CRAN.R-project.org/package=gamlss). These data contain 86 annual maximum seas levels recorded at Fremantle, Australia during 1987-1989. In addition to the year of each sea level, we have available the value of the Southern Oscillation Index (SOI). We use the `fitGEV()` function provided in `gamlssx` to fit a model to these data that is similar to the first one fitted, to the same data, in Gavin Simpson's blog post. 

The `fitGEV()` function calls the function `gamlss::gamlss()`, which offers 3 fitting algorithms: `RS` (Rigby and Stasinopoulos), `CG` (Cole and Green) and `mixed` (`RS` initially followed by `CG`). In the code below, we use the default `RS` algorithm. `fitGEV()` offers 2 scoring methods to calculate the weights used in the algorithm. Here, we use the default, Fisher's scoring, based on the expected Fisher information. The code below does not do justice to the functionality of the `gamlss` package.  See the [GAMLSS books](https://www.gamlss.com/books-vignettes/) for more information. 

```{r, message = FALSE}
# Load gamlss, for the function term.plot()
library(gamlss)
# Load gamlssx
library(gamlssx)
# Transform Year so that it is centred on 0 
fremantle <- transform(fremantle, cYear = Year - median(Year))
```

```{r, fig.alt = "Plot of sea level against year"}
# Plot sea level against year and against SOI
plot(fremantle$Year, fremantle$SeaLevel, xlab = "year", ylab = "sea level (m)")
```

```{r, fig.alt = "Plot of sea level against SOI"}
plot(fremantle$SOI, fremantle$SeaLevel, xlab = "SOI", ylab = "sea level (m)")
```

```{r}
# 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)
```

```{r, fig.alt = "Residual diagnostic plots"}
# Model diagnostic plots
plot(mod)
```

```{r, fig.alt = "Term plots for parameter mu for cYear and SOI"}
# 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)
```

```{r, fig.alt = "Term plots for parameter sigma for cYear and SOI"}
# Scale sigma
term.plot(mod, what = "sigma", rug = TRUE, pages = 1)
```

## Installation

To get the current released version from CRAN:

```{r installation, eval = FALSE}
install.packages("gamlssx")
```

Owner

GitHub Events

Issues and Pull Requests

Last synced: 6 months ago

Packages