tests/testthat/test-owtemps-CB2007.R

In paulnorthrop/oola: Objected-oriented Loglikelihood Adjustment

context("owtemp_CB2007")

# Check that oola (using evd::fgev) and chandwich agree on the analyses of the
# Oxford-Worthing temperature data in Chandler and Bate (2007).

my_tol <- 1e-4

# -------------------------- GEV model, owtemps data -----------------------
# ------------ following Section 5.2 of Chandler and Bate (2007) -----------

# Contributions to the independence loglikelihood
gev_loglik <- function(pars, data) {
  o_pars <- pars[c(1, 3, 5)] + pars[c(2, 4, 6)]
  w_pars <- pars[c(1, 3, 5)] - pars[c(2, 4, 6)]
  if (o_pars[2] <= 0 | w_pars[2] <= 0) return(-Inf)
  o_data <- data[, "Oxford"]
  w_data <- data[, "Worthing"]
  check <- 1 + o_pars[3] * (o_data - o_pars[1]) / o_pars[2]
  if (any(check <= 0)) return(-Inf)
  check <- 1 + w_pars[3] * (w_data - w_pars[1]) / w_pars[2]
  if (any(check <= 0)) return(-Inf)
  o_loglik <- chandwich::log_gev(o_data, o_pars[1], o_pars[2], o_pars[3])
  w_loglik <- chandwich::log_gev(w_data, w_pars[1], w_pars[2], w_pars[3])
  return(o_loglik + w_loglik)
}

test_data <- chandwich::owtemps

got_evd <- requireNamespace("evd", quietly = TRUE)

if (got_evd) {
  # Fit all the models using evd::fgev
  y <- c(test_data[, "Oxford"], test_data[, "Worthing"])
  year <- rep(1:(length(y) / 2), 2)
  x <- rep(c(1, -1), each = length(y) / 2)
  owfit_s <- evd::fgev(y, nsloc = x)
  init_s <- c(owfit_s$estimate[1:3], 0, owfit_s$estimate[4], 0)
  owfit_t <- evd::fgev(y)
  init_t <- c(owfit_t$estimate[1], 0, owfit_t$estimate[2], 0,
              owfit_t$estimate[3], 0)
  par_names <- c("mu[0]", "mu[1]", "sigma[0]", "sigma[1]", "xi[0]", "xi[1]")
  small <- chandwich::adjust_loglik(gev_loglik, data = test_data,
                                    init = init_s,
                                    par_names = par_names,
                                    fixed_pars = c("sigma[1]", "xi[1]"))
  tiny <- chandwich::adjust_loglik(gev_loglik, data = test_data,
                                   init = init_t,
                                   par_names = par_names,
                                   fixed_pars = c("mu[1]", "sigma[1]", "xi[1]"))

  # Note: we would need to adjust nobs in small and tiny (it's half of that in
  # oola_small and oola_tiny)
  attr(small, "nobs") <- attr(small, "nobs") * 2
  attr(tiny, "nobs") <- attr(tiny, "nobs") * 2

  # Use the MLE from owfit_s to aid the comparison
  # Set use_vcov = FALSE so that we estimate H inside adjust_loglik()
  oola_small <- alogLik(owfit_s, cluster = year, cadjust = FALSE,
                        use_vcov = FALSE)
  test_that("owtemps, small, logLiks: oola vs. chandwich", {
    testthat::expect_equivalent(logLik(oola_small), logLik(small),
                                tolerance = my_tol)
  })
  test_that("owtemps, small, adjSEs: oola vs. chandwich", {
    testthat::expect_equivalent(attr(oola_small, "adjSE"),
                                attr(small, "adjSE"), tolerance = my_tol)
  })
  test_that("owtemps, small, adjVC: oola vs. chandwich", {
    testthat::expect_equivalent(attr(oola_small, "adjVC"),
                                attr(small, "adjVC"), tolerance = my_tol)
  })

  # Use the MLE from owfit_t to aid the comparison
  # Set use_vcov = FALSE so that we estimate H inside adjust_loglik()
  oola_tiny <- alogLik(owfit_t, cluster = year, cadjust = FALSE,
                       use_vcov = FALSE)
  test_that("owtemps, tiny, logLiks: oola vs. chandwich", {
    testthat::expect_equivalent(logLik(oola_tiny), logLik(tiny),
                                tolerance = my_tol)
  })
  test_that("owtemps, tiny, adjSEs: oola vs. chandwich", {
    testthat::expect_equivalent(attr(oola_tiny, "adjSE"),
                                attr(tiny, "adjSE"), tolerance = my_tol)
  })
  test_that("owtemps, tiny, adjVC: oola vs. chandwich", {
    testthat::expect_equivalent(attr(oola_tiny, "adjVC"),
                                attr(tiny, "adjVC"), tolerance = my_tol)
  })

  # Test anova methods
  res1 <- anova(small, tiny)
  res2 <- anova(oola_small, oola_tiny)
  test_that("owtemps, anova, small vs tiny", {
    testthat::expect_equivalent(res1, res2, tolerance = my_tol)
  })
}