tests/testthat/test-logLik.R

In paulnorthrop/oola: Objected-oriented Loglikelihood Adjustment

context("logLik")

# Check that the same loglikelihood values and df are returned from
# 1. the standard logLik S3 method
# 2. my logLikVec S3 method, after summing to evaluate the log-likelihood
#    using logLik.logLikVec()
# 3. logLik.chandwich()

# Note: logLik.chandwich() in version 1.1.0 of chandwich on CRAN didn't return
#       nobs as an attribute, which it should.  This has been fixed, but, at
#       least for the moment we use expect_equivalent(), rather than
#       expect_equal(), to avoid the problem that the length of
#       the attributes return from logLik.chandwich is 2 rather than 3.

# Binomial GLM, vector response

if (requireNamespace("AER", quietly = TRUE)) {
  data("Affairs", package = "AER")
  fm_probit <- stats::glm(I(affairs > 0) ~ age + yearsmarried + religiousness +
                          occupation + rating, data = Affairs,
                          family = binomial(link = "probit"))
  adj_binom <- alogLik(fm_probit)

  test_that("binom GLM, vector response: logLik.glm() vs. logLik.chandwich()", {
    testthat::expect_equivalent(logLik(fm_probit), logLik(adj_binom))
  })
  test_that("binom GLM, vector response: logLik.glm() vs. logLik.binomglm()", {
    testthat::expect_equal(logLik(fm_probit), logLik(logLikVec(fm_probit)))
  })
}

# Binomial GLM, matrix response

lrfit <- stats::glm(cbind(using, notUsing) ~ age + education + wantsMore,
                    family = binomial, data = cuse)
adj_binom2 <- alogLik(lrfit)

test_that("binom GLM, matrix response: logLik.glm() vs. loglik.chandwich()", {
  testthat::expect_equivalent(logLik(lrfit), logLik(adj_binom2))
})
test_that("binom GLM, matrix response: logLik.glm() vs. loglik.binomglm()", {
  testthat::expect_equal(logLik(lrfit), logLik(logLikVec(lrfit)))
})

# Poisson GLM

counts <- c(18,17,15,20,10,20,25,13,12)
outcome <- gl(3,1,9)
treatment <- gl(3,3)
d.AD <- data.frame(treatment, outcome, counts)
glm.D93 <- stats::glm(counts ~ outcome + treatment, family = poisson,
                      data = d.AD)
adj_pois <- alogLik(glm.D93)

test_that("pois GLM: logLik.glm() vs. loglik.chandwich()", {
  testthat::expect_equivalent(logLik(glm.D93), logLik(adj_pois))
})
test_that("pois GLM: logLik.glm() vs. loglik.poisglm()", {
  testthat::expect_equal(logLik(glm.D93), logLik(logLikVec(glm.D93)))
})

# evd::fgev()

if (requireNamespace("evd", quietly = TRUE)) {
  uvdata <- evd::rgev(100, loc = 0.13, scale = 1.1, shape = 0.2)
  trend <- (-49:50)/100

  # Estimate trend
  M1 <- evd::fgev(uvdata, nsloc = trend)
  adj_fgev <- alogLik(M1)
  temp <- M1
  class(temp) <- "gev"

  test_that("evd::fgev: logLik.gev() vs. loglik.chandwich()", {
    # Note: evd::logLik.evd returns non-standard attributes (no nobs)
    # Therefore, use expect_equivalent(), rather than expect_equal()
    testthat::expect_equivalent(logLik(M1), logLik(adj_fgev))
  })
  test_that("evd::fgev: logLik.chandwich() vs. loglikVec.gev()", {
    testthat::expect_equivalent(logLik(adj_fgev), logLik(logLikVec(temp)))
  })

  # Set shape = 0
  M2 <- evd::fgev(uvdata, nsloc = (-49:50)/100, shape = 0)
  adj_fgev <- alogLik(M2)
  temp <- M2
  class(temp) <- "gev"

  test_that("evd::fgev: logLik.gev() vs. loglik.chandwich()", {
    # Note: evd::logLik.evd returns non-standard attributes (no nobs)
    # Therefore, use expect_equivalent(), rather than expect_equal()
    testthat::expect_equivalent(logLik(M2), logLik(adj_fgev))
  })
  test_that("evd::fgev: logLik.chandwich() vs. loglikVec.gev()", {
    testthat::expect_equivalent(logLik(adj_fgev), logLik(logLikVec(temp)))
  })

}