evgam (version 1.0.2)

custom.family.evgam: Custom distributions with evgam

Description

Users can supply code to fit an arbitrary distribution with evgam using family = "custom". See Details for more information and the example below, which demonstrates fitting of the Gumbel distribution.

Details

Users should supply a list to evgam called custom.fns, which which comprises the following functions: d0, which evaluates the negative log-likelihood; d120, which evaluates its first and second derivatives; and, optionally, d340, which evaluates its third and fourth derivatives. The list may also contain q, which evaluates the custom distribution's quantile function, for use with predict, prob = ...). The list may also contain unlink, which gives functions to reverse any link functions (unlink) used with linear predictors, so that predict(..., type = "response") works meaningfully, and unlink may also have an attribute deriv, which gives derivatives of unlink functions, which allows predict(..., type = "response", std.err = TRUE) to return meaningful values.

Trying to mimic the Gumbel example below for your chosen distribution is almost certainly the easiest was to get family = "custom" to do what you want.

References

Youngman, B. D. (2022). evgam: An R Package for Generalized Additive Extreme Value Modules. Journal of Statistical Software. To appear. tools:::Rd_expr_doi("10.18637/jss.v103.i03")

See Also

family.evgam

Examples

# Gumbel custom likelihood

# negative log likelihood
gumd0 <- function(pars, likdata) {

# this is needed, and should be fine
pars <- split(pars, likdata$idpars)

# likdata$X should be set up by evgam
mu <- likdata$X[[1]] %*% pars[[1]]
lsigma <- likdata$X[[2]] %*% pars[[2]]
y <- likdata$y

y <- (y - mu) * exp(-lsigma)
nllh <- sum(lsigma - y + exp(y))

return(nllh)

}

# first and second derivatives of neg. log lik. in order
# d_mu, d_lsigma, d_{mu, mu}, d_{mu, lsigma}, d_{lsigma, lsigma}
gumd12 <- function(pars, likdata) {

# this is needed, and should be fine
pars <- split(pars, likdata$idpars)

# likdata$X should be set up by evgam
mu <- likdata$X[[1]] %*% pars[[1]]
lsigma <- likdata$X[[2]] %*% pars[[2]]
y <- likdata$y

out <- matrix(0, length(y), 5)

ee1 <- exp(lsigma)
ee2 <- y - mu
ee3 <- ee2/ee1
ee4 <- exp(ee3)
ee5 <- (ee3 + 1) * ee4
ee6 <- 1 - ee4

# first derivatives
out[, 1] <- ee6/ee1
out[, 2] <- ee6 * ee2/ee1 + 1

# second derivatives
out[, 3] <- ee4/ee1^2
out[, 4] <- -((1 - ee5)/ee1)
out[, 5] <- (ee5 - 1) * ee2/ee1

return(out)

}

gum_fns <- list(d0 = gumd0, d120 = gumd12, d340 = NULL)

unlink <- list(NULL, function(x) exp(x))
attr(unlink[[2]], "deriv") <- unlink[[2]]
qgumbel <- function(p, pars1, pars2) pars1 - exp(pars2) * log(-log(p))

gum_fns$q <- qgumbel
gum_fns$unlink <- unlink


# \donttest{

data(COprcp)
COprcp <- cbind(COprcp, COprcp_meta[COprcp$meta_row,])
COprcp$year <- format(COprcp$date, "%Y")
COprcp_gev <- aggregate(prcp ~ year + meta_row, COprcp, max)
COprcp_gev <- cbind(COprcp_gev, COprcp_meta[COprcp_gev$meta_row,])

inits <- sqrt(6 * var(COprcp_gev$prcp) / pi^2)
inits <- c(mean(COprcp_gev$prcp) - inits[1] * 0.5772156649, log(inits[1]))
fmla_gum <- list(location = prcp ~ s(lon, lat) + s(elev), logscale = ~ s(lon, lat))
m <- evgam(fmla_gum, data = COprcp_gev, family = 'custom', custom.fns = gum_fns, 
trace = 2, inits = inits)

predict(m, COprcp_gev[1:10,])
predict(m, COprcp_gev[1:10,], type = 'response')
predict(m, COprcp_gev[1:10,], prob = .99)

# }