In paulnorthrop/oola: Objected-oriented Loglikelihood Adjustment

knitr::opts_chunk$set(comment = "#>", collapse = TRUE)

For an overview of the oola package see the vignette oola: Object-oriented loglikelihood adjustment.

The evd package @evd

The function fgev in the evd package fits a Generalized Extreme Value (GEV) distribution to data using maximum likelihood estimation. A linear regression effect of covariates on the location parameter may be included.

@CB2007 @sandwich

Unless stated otherwise, the examples below are taken from the help files of the relevant function.

NOTE: add illustrations for what oola gives re loglikelihoods.

library(oola)
got_evd <- requireNamespace("evd", quietly = TRUE)
if (got_evd) {
  library(evd)
  uvdata <- evd::rgev(100, loc = 0.13, scale = 1.1, shape = 0.2)
  M1 <- evd::fgev(uvdata, nsloc = (-49:50)/100)
}  

Setting logLikVec

We create a function to calculate the

logLikVec.gev <- function(object, pars = NULL, ...) {
  if (!missing(...)) {
  warning("extra arguments discarded")
  }
  # If the parameter estimates have not been provided in pars then extract
  # them from the fitted object
  if (is.null(pars)) {
    pars <- coef(object)
  }
  n_pars <- length(pars)
  mu <- pars[1]
  sigma <- pars[n_pars - 1]
  xi <- pars[n_pars]
  if (n_pars > 3) {
    mu_reg <- pars[2:(n_pars - 2)]
    mu <- mu + as.matrix(object$nsloc) %*% mu_reg
  }
  # Calculate the weighted loglikelihood contributions
  if (sigma <= 0) {
    val <- -Inf
  } else {
    val <- evd::dgev(object$data, loc = as.vector(mu), scale = sigma,
                     shape = xi, log = TRUE)
  }
  # Return the usual attributes for a "logLik" object
  attr(val, "nobs") <- length(object$data)
  attr(val, "df") <- n_pars
  class(val) <- "logLikvec"
  return(val)
}

The evd package has a logLik method for objects inheriting from class "evd". oola has a loglik method for objects of class logLikVec, which simply sums the loglikelihood contributions and retains the attributes of the object. We check that the values of the maximized loglikelihoods returned from these two approaches agree.

class(M1)
logLik(M1)
logLik(logLikVec(M1))

Setting nobs, coef and vcov (if required)

The object whose loglikelihood is to be adjusted must have S3 methods nobs, coef and vcov methods available. If they are not available then the user must create them. In the current case, the evd package provides a vcov method for objects returned from evd::fgev() but does not provide nobs or coef methods. We provide these below and, for completeness, also provide a vcov method.

nobs.gev <- function(object, ...) {
  return(object$n)

coef.gev <- function(object, ...) {
  return(object$estimate)
}

vcov.gev <- function(object, ...) {
  return(object$var.cov)
}

Adjusting the log-likelihood

if (got_evd) {
  adj_fgev <- alogLik(M1)
  summary(adj_fgev)
}

References

paulnorthrop/oola documentation built on May 12, 2019, 10:52 a.m.