tests/testthat/test-hskewL.R

In gnm: Generalized Nonlinear Models

# some tests from tests/hskewL.R in logmult package

test_that("gnmFit handles extra linear constraints with eliminate", {
    if (requireNamespace("logmult")){
        data(ocg1973, package = "logmult")
        
        tab <- array(ocg1973, dim=c(nrow(ocg1973), ncol(ocg1973), 2))
        
        # model contains extra linear constraints
        model <- logmult::hmskewL(tab[5:1, 5:1,], weighting="uniform", start=NA)
        
        ass <- model$assoc
        
        # First score for Farmers is slightly different from original article
        expect_equal(round(ass$row[,,1] * sqrt(ass$phi[1,1]), d=2)[5:1,],
                     matrix(c(-0.08, -0.2, -0.23, -0.11, 0.61,
                              0.34,  0.3, -0.13, -0.51, 0), 5, 2),
                     ignore_attr = TRUE)
        expect_equal(round(ass$row[,,1] * sqrt(ass$phi[2,1]), d=2)[5:1,],
                     matrix(c(-0.08, -0.2, -0.23, -0.11, 0.61,
                              0.34,  0.3, -0.13, -0.51, 0), 5, 2),
                     ignore_attr = TRUE)

    }
})

test_that("gnmFit handles computationally singular initial linear fit ", {
    if (requireNamespace("logmult")){
        # EGP class of cohabiting spouses where one is 30-60 
        # (last occupation for inactive persons)
        # French Labour Force Surveys, 1969 and 2011
        tab2 <- structure(
            c(261L, 43L, 21L, 5L, 7L, 26L, 5L, 16L, 17L, 7L, 1L,
              483L, 394L, 215L, 53L, 58L, 117L, 37L, 185L, 232L, 104L, 11L,
              565L, 457L, 528L, 139L, 116L, 201L, 19L, 469L, 788L, 368L, 14L,
              148L, 195L, 399L, 306L, 96L, 213L, 50L, 321L, 1327L, 1344L, 123L,
              17L, 9L, 11L, 3L, 33L, 37L, 4L, 12L, 13L, 11L, 1L, 165L, 70L,
              105L, 77L, 573L, 878L, 55L, 78L, 188L, 181L, 11L, 9L, 3L, 34L,
              7L, 23L, 36L, 1918L, 13L, 88L, 228L, 87L, 26L, 10L, 16L, 3L,
              0L, 1L, 2L, 27L, 31L, 24L, 1L, 88L, 115L, 191L, 57L, 50L, 98L,
              32L, 201L, 552L, 493L, 19L, 49L, 115L, 317L, 122L, 58L, 110L,
              38L, 316L, 1301L, 1622L, 66L, 0L, 3L, 9L, 6L, 4L, 13L, 15L, 7L,
              56L, 135L, 143L, 919L, 189L, 54L, 32L, 74L, 64L, 19L, 113L, 86L,
              40L, 3L, 875L, 519L, 183L, 97L, 129L, 129L, 62L, 329L, 343L,
              195L, 12L, 513L, 330L, 271L, 126L, 188L, 145L, 70L, 382L, 578L,
              388L, 15L, 250L, 236L, 180L, 217L, 126L, 155L, 52L, 356L, 965L,
              634L, 42L, 28L, 8L, 10L, 5L, 59L, 14L, 2L, 20L, 30L, 6L, 1L,
              61L, 30L, 20L, 7L, 52L, 106L, 4L, 23L, 38L, 17L, 2L, 4L, 1L,
              2L, 4L, 6L, 3L, 135L, 4L, 7L, 8L, 2L, 53L, 19L, 7L, 3L, 8L, 6L,
              3L, 39L, 31L, 22L, 2L, 28L, 44L, 34L, 20L, 19L, 16L, 11L, 64L,
              165L, 105L, 9L, 41L, 35L, 42L, 52L, 20L, 38L, 10L, 111L, 299L,
              309L, 12L, 1L, 3L, 2L, 2L, 0L, 4L, 25L, 4L, 32L, 19L, 16L),
            .Dim = c(11L, 11L, 2L), class="table",
            .Dimnames = structure(
                list(H = c("I", "II", "IIIa", "IIIb", "IVa",
                           "IVb", "IVc", "V", "VI", "VIIa", "VIIb"),
                     F = c("I", "II", "IIIa", "IIIb", "IVa", "IVb",
                           "IVc", "V", "VI", "VIIa", "VIIb"),
                     T = c("1969", "2011")),
                .Names = c("M", "W", "T")))
        
        # requires use of ridge in glm.fit.e 
        model2 <- logmult::hmskewL(tab2, start=NA)
        
        expect_equal(round(c(model2$assoc$phi), 2), c(0.18, 0.04, 0.18, 0.04))
        expect_equal(round(c(model2$assoc$row), 2),
                     c(1.97, 1.38,  0.05, -0.24, -1.02, 0.57, -0.79, -0.16,
                       -0.14, -1.32, -1.56,
                       0,   -0.03, -0.94,  0.98, -0.10, 0.71,  2.65, -0.86,
                       -0.66, -0.77,  0.79))
    }
})