tests/testthat/test-glm.R

context("enrichment of glm objects")

library("numDeriv")
library("MASS")

## A Gamma example, from McCullagh & Nelder (1989, pp. 300-2)
clotting <- data.frame(
    u = c(5,10,15,20,30,40,60,80,100, 5,10,15,20,30,40,60,80,100),
    conc = c(118,58,42,35,27,25,21,19,18,69,35,26,21,18,16,13,12,12),
    lot = factor(c(rep(1, 9), rep(2, 9))))
mod1 <- glm(conc ~ lot*log(u), data = clotting, family = Gamma)

tol <- 1e-05
test_that("implementation of the scores corresponds to that of the observed information [glm]",
{
    enriched_mod1 <- enrich(mod1, with = "auxiliary functions")
    disp <- summary(enriched_mod1)$disp
    info_appr <- -jacobian(function(coefs) {
        p <- length(coefs)
        enriched_mod1$auxiliary_functions$score(coefs[1:(p-1)],
                                          dispersion = coefs[p])
    }, c(coef(enriched_mod1), disp))
    info_exac <- enriched_mod1$auxiliary_functions$information(coef(enriched_mod1),
                                                         dispersion = disp,
                                                         type = "observed")
    expect_equal(sum(abs(info_appr - info_exac)), 0, tolerance = 0.01)
})

test_that("ML estimate of gamma dispersion from enrichwith is numerically the same to that from MASS::gamma.dispersion", {
    enriched_mod1 <- enrich(mod1, with = "mle of dispersion")
    expect_equal(unname(enriched_mod1$dispersion_mle),
                 gamma.dispersion(mod1),
                 tolerance = tol)
})

test_that("standard errors from expected information matrix match those from summary", {
    dispersion_pearson <- summary(mod1)$dispersion
    enriched_mod1 <- enrich(mod1, with = "expected information", dispersion = dispersion_pearson)
    stderrs1 <- sqrt(diag(solve(enriched_mod1$expected_information)))
    stderrs2 <- coef(summary(mod1))[, 2]
    expect_equal(stderrs1[-length(stderrs1)], stderrs2)
})


test_that("enrichwith handles aliasing correctly [glm]",
{

    umod1 <- update(mod1, . ~ . + I(2 * log(u)) + I(3 * log(u)))
    na_coefs <- is.na(coef(umod1))
    scorefun <- get_score_function(umod1)
    infofun <- get_information_function(umod1)
    biasfunA <- get_bias_function(umod1)
    biasfunB <- get_bias_function(mod1)
    expect_true(all(is.na(scorefun()[na_coefs])))
    expect_true(all(is.na(infofun()[na_coefs, ])))
    expect_true(all(is.na(infofun()[, na_coefs])))
    expect_true(all(biasfunA()[names(biasfunB())] == biasfunB()))

})

## A Poisson example
## Dobson (1990) Page 93: Randomized Controlled Trial :
counts <- c(18,17,15,20,10,20,25,13,12)
outcome <- gl(3,1,9)
treatment <- gl(3,3)
mod2 <- glm(counts ~ outcome + treatment, family = poisson())

test_that("expected information matrix from enrich with is equal to that coming from the qr decomposition component of the object [glm]", {
    enriched_mod2 <- enrich(mod2, with = "auxiliary functions")
    expect_equal(enriched_mod2$auxiliary_functions$information(coef(mod2), type = "expected"),
                 crossprod(qr.R(mod2$qr)),
                 tolerance = 0.0001,
                 check.attributes = FALSE)
})
ikosmidis/enrichwith documentation built on Jan. 1, 2020, 9:44 a.m.