tests/testthat/test-dropouts.R

In OrthoPanels: Dynamic Panel Models with Orthogonal Reparameterization of Fixed Effects

context('dropouts')

## Generate test data without littering the environment with temporary
## variables
x <- NULL
y <- NULL
d <- NULL
o <- NULL
local({
    suppressWarnings(RNGversion("3.5.0"))
    set.seed(123, kind = "Mersenne-Twister", normal.kind = "Inversion")

    N <- 10
    T <- 3
    
    dd <- generate_data(N=N, T=T, drop_outs=2)
    x <<- dd$x
    y <<- dd$y
    d <<- data.frame(i = rep(seq_len(N), each=T+1),
                     t = rep(seq_len(T+1), N),
                     x = c(x),
                     y = c(y))

    suppressWarnings(RNGversion("3.5.0"))
    set.seed(123, kind = "Mersenne-Twister", normal.kind = "Inversion")
    o <<- opm(y~x, d, n.samp = 10)
})

## Sanity check
test_that('data', {
    expect_equal(x, array(c(-1.31047564655221, 0.474081797439462, -1.81782370598685,
-0.323535778523186, -0.23017748948328, 0.359813827057364, -0.217974914658295,
-0.295071482992271, 2.30870831414912, 1.15077145059405, -0.27600444830724,
1.64512566104502, -0.679491608575424, -0.63931728405488, NA,
0.128133487533042, 0.129287735160946, -0.555841134754075, -0.625039267849257,
0.821581081637487, 2.46506498688328, 2.53691313680308, NA, 1.43864025410009,
-0.289083794010798, -0.252149521770761, 0.0877870444945247, -0.196082346462411,
-1.26506123460653, -1.96661715662964, 0.153373117836515, -0.0619117105767217,
0.0631471481064739, 1.45135590156369, -0.388136937011948, 0.444037336260083,
-1.19566197009996, -1.22279140772793, 0.503814921069927, -1.13047100101238
),
                           dim = c(4, 1, 10)))
    expect_equal(y, matrix(c(-2.34994480219662, -1.68461298838382, -2.37157886442545,
-2.83858848753085, -0.323006022761239, -0.0101428532006405, -0.61638233703877,
-2.76489578565633, 0.888957805506297, 1.97699417075886, 1.51728747755639,
3.58694509376296, 0.829210161050801, 0.463548722512418, NA, -2.9082573370641,
1.27260586588546, 0.132611379906427, -1.31800517044699, -0.936220660872111,
1.10942391023829, 4.33963912795022, NA, 4.24739227175083, -1.54742673230447,
-3.44854093126784, -2.23216716475723, -2.49889776266083, -1.09918597092649,
-0.948287814141993, -0.344453121422235, -1.42390012825401, 1.81153869238955,
2.75530154082123, 3.10584976978438, 2.95624703277138, -1.68120005152181,
-2.2360541608809, 0.183965065721659, -1.61214433008441),
                           nrow = 4, ncol = 10))
})





test_that('default', {
    suppressWarnings(RNGversion("3.5.0"))
    set.seed(123, kind = "Mersenne-Twister", normal.kind = "Inversion")
    expect_equal(opm(x, y, n = 10),
                 structure(list(samples = list(rho = c(0.847, 0.466,
                                                       0.767, 0.763, 0.577, 0.973, 0.688, 0.553, 0.672, 0.735),
                                               sig2 = c(0.577161939376618,
                                                        1.37678799790679, 0.547297492131444, 0.434998550554306, 0.771225976041554,
                                                        0.723112120326027, 0.920962982724747, 1.3047157196198, 1.41232132729643,
                                                        2.00647541468943),
                                               beta = matrix(c(0.902648414988294, 0.0632238525152027,
                                                               0.873006012775479, 0.868237559365092, 0.423891867329579, 1.05821070795326,
                                                               0.605189520563423, 0.564474185116966, 0.516742369586711, 0.289816166868861
                                                               ), 10, 1)),
                                fitted.values = matrix(c(0.449620185553199,
                                                         -0.417375571022443, -0.032244614530756, 0.235680717369341, 0.120373602143934,
                                                         -0.356054319513275, -0.225294843952018, -0.285580962616997, 0.510875806569015,
                                                         0, 0, 0, 0.629744904853414, -0.221830444311809, -0.407914460541605,
                                                         0, 0, 0, 0.53607998944547, -0.617757031219903, 0.0816770417744323,
                                                         -0.99937557210787, 0.347826306663529, 0.651549265444341, 0.024172129930576,
                                                         -0.38013460848074, 0.355962478550163, -0.665350858340203, -0.0503551787527371,
                                                         0.71570603709294), 3, 10),
                                residuals = matrix(c(0.164026939509687,
                                                     0.3440568200437, -0.508083759553388, 0.884650088061933, 0.393717719449211,
                                                     -1.27836780751114, -0.158119899315194, -0.557540473852684, 0.715660373167878,
                                                     0, 0, 0, 0.210071292190572, -0.388969908997625, 0.178898616807053,
                                                     0, 0, 0, -1.25808563448468, 1.11212515269131, 0.145960481793371,
                                                     0.956634779238624, 0.213267593186983, -1.16990237242561, -0.208003370235008,
                                                     0.546851597139455, -0.338848226904447, -0.349292160792813, 1.45573138622228,
                                                     -1.10643922542946), 3, 10),
                                df.residual = sum(c(3, 3, 3, 1, 3, 1, 3, 3, 3, 3)) - 2,
                                logLik = -1.12486246022849,
                                design = "unbalanced (with dropouts)",
                                call = quote(opm(x = x, y = y, n.samp = 10)),
                                .Environment = environment(),
                                time.indicators = FALSE),
                           class = "opm"))
})


test_that('late-arrivals', {
    ## shift 4th and 6th subject so that they start late
    ## this shouldn't change the result
    x[,,4] <- c(NA, x[1:3, , 4])
    y[,4] <- c(NA, y[1:3, 4])
    x[,,6] <- c(NA, x[1:3, , 6])
    y[,6] <- c(NA, y[1:3, 6])
    
    suppressWarnings(RNGversion("3.5.0"))
    set.seed(123, kind = "Mersenne-Twister", normal.kind = "Inversion")
    expect_equal(opm(x, y, n = 10),
                 structure(list(samples = list(rho = c(0.847, 0.466,
                                                       0.767, 0.763, 0.577, 0.973, 0.688, 0.553, 0.672, 0.735),
                                               sig2 = c(0.577161939376618,
                                                        1.37678799790679, 0.547297492131444, 0.434998550554306, 0.771225976041554,
                                                        0.723112120326027, 0.920962982724747, 1.3047157196198, 1.41232132729643,
                                                        2.00647541468943),
                                               beta = matrix(c(0.902648414988294, 0.0632238525152027,
                                                               0.873006012775479, 0.868237559365092, 0.423891867329579, 1.05821070795326,
                                                               0.605189520563423, 0.564474185116966, 0.516742369586711, 0.289816166868861
                                                               ), 10, 1)),
                                fitted.values = matrix(c(0.449620185553199,
                                                         -0.417375571022443, -0.032244614530756, 0.235680717369341, 0.120373602143934,
                                                         -0.356054319513275, -0.225294843952018, -0.285580962616997, 0.510875806569015,
                                                         0, 0, 0, 0.629744904853414, -0.221830444311809, -0.407914460541605,
                                                         0, 0, 0, 0.53607998944547, -0.617757031219903, 0.0816770417744323,
                                                         -0.99937557210787, 0.347826306663529, 0.651549265444341, 0.024172129930576,
                                                         -0.38013460848074, 0.355962478550163, -0.665350858340203, -0.0503551787527371,
                                                         0.71570603709294), 3, 10),
                                residuals = matrix(c(0.164026939509687,
                                                     0.3440568200437, -0.508083759553388, 0.884650088061933, 0.393717719449211,
                                                     -1.27836780751114, -0.158119899315194, -0.557540473852684, 0.715660373167878,
                                                     0, 0, 0, 0.210071292190572, -0.388969908997625, 0.178898616807053,
                                                     0, 0, 0, -1.25808563448468, 1.11212515269131, 0.145960481793371,
                                                     0.956634779238624, 0.213267593186983, -1.16990237242561, -0.208003370235008,
                                                     0.546851597139455, -0.338848226904447, -0.349292160792813, 1.45573138622228,
                                                     -1.10643922542946), 3, 10),
                                df.residual = sum(c(3, 3, 3, 1, 3, 1, 3, 3, 3, 3)) - 2,
                                logLik = -1.12486246022849,
                                design = "unbalanced (with dropouts)",
                                call = quote(opm(x = x, y = y, n.samp = 10)),
                                .Environment = environment(),
                                time.indicators = FALSE),
                           class = "opm"))
    
    ## shift 4th and 6th subject so that they start late, and don't
    ## have an entry for t=1
    ## this shouldn't change the result
    d[d$i == 4,'t'] <- 2:5
    d[d$i == 6,'t'] <- 2:5
    d <- subset(d, t != 5)
    o$.Environment = environment()
    suppressWarnings(RNGversion("3.5.0"))
    set.seed(123, kind = "Mersenne-Twister", normal.kind = "Inversion")
    expect_equal(opm(y ~ x, d, n.samp = 10), o)
})


test_that('formula', {
    suppressWarnings(RNGversion("3.5.0"))
    set.seed(123, kind = "Mersenne-Twister", normal.kind = "Inversion")
    expected <- o
    expected$.Environment <- environment()
    expect_equal(opm(y~x, d, n.samp = 10), expected)
    
    ## different order for i and t columns
    d$i <- letters[d$i]
    d$t <- d$t*1000
    d <- d[sample(nrow(d)), ]
    dimnames(expected$residuals) <- list(t = 2:4 * 1000,
                                         i = letters[1:10])
    set.seed(123)
    expect_equal(opm(y~x, d, n.samp = 10), expected)
})

test_that('DIC', {
    expect_equal(DIC(o), 15.1580003661899)
})