1library("lme4")
2library("testthat")
3testLevel <- if (nzchar(s <- Sys.getenv("LME4_TEST_LEVEL"))) as.numeric(s) else 1
4
5set.seed(101)
6dd <- expand.grid(f1 = factor(1:3),
7                  f2 = LETTERS[1:2], g=1:9, rep=1:15,
8          KEEP.OUT.ATTRS=FALSE)
9mu <- 5*(-4 + with(dd, as.integer(f1) + 4*as.numeric(f2)))
10dd$y <- rnbinom(nrow(dd), mu = mu, size = 0.5)
11
12## mimic glmer.nb protocol
13if (testLevel>1) {
14test_that("most messages suppressed", {
15    expect_message(glmer.nb(y ~ f1 + (1|g), data=dd[1:10,]),
16                   "singular")
17})
18
19test_that("ok with negative.binomial masking", {
20    negative.binomial <- function() {}
21    ## use shortened version of data for speed ...
22    m.base <- glmer.nb(y ~ f1 + (1|g), data=dd[1:200,])
23    expect_is(m.base,"merMod")
24})
25
26test_that("ok with Poisson masking", {
27    poisson <- NA
28    ## use shortened version of data for speed ...
29    m.base <- glmer.nb(y ~ f1 + (1|g), data=dd[1:200,])
30    expect_is(m.base,"merMod")
31    rm(poisson)
32})
33
34if (testLevel>2) {
35context("testing glmer refit")
36
37test_that("glmer refit", {
38            ## basic Poisson fit
39            m.base <- glmer(y ~ f1*f2 + (1|g), data=dd, family=poisson)
40            expect_equal(m.base@beta,(m.base.r <- refit(m.base))@beta,
41                         tolerance = 1e-5)
42
43            th <- lme4:::est_theta(m.base,limit=20,eps=1e-4,trace=FALSE)
44            th0 <- structure(0.482681268108477, SE = 0.0244825021248148)
45            th1 <- structure(0.482681277470945)
46            th2 <- 0.482681268108477
47            th3 <- 0.4826813
48            ## NB update with raw number
49            m.numth1 <- update(m.base,family=MASS::negative.binomial(theta=0.4826813))
50            expect_equal(m.numth1@beta,(m.numth1.r <- refit(m.numth1))@beta)
51
52            ## strip NB value
53            m.symth4 <- update(m.base,family=MASS::negative.binomial(theta=c(th)))
54            expect_equal(m.symth4@beta,(m.symth4.r <- refit(m.symth4))@beta)
55
56            ## IDENTICAL numeric value to case #1 above
57            m.symth6 <- update(m.base,family=MASS::negative.binomial(theta=th3))
58            expect_equal(m.symth6@beta,(m.symth6.r <- refit(m.symth6))@beta)
59
60            ## standard NB update with computed theta from est_theta (incl SE attribute)
61            m.symth <- update(m.base,family=MASS::negative.binomial(theta=th))
62            expect_equal(m.symth@beta,(m.symth.r <- refit(m.symth))@beta)
63
64            ## NB update with equivalent value
65            m.symth2 <- update(m.base,family=MASS::negative.binomial(theta=th0))
66            expect_equal(m.symth2@beta,(m.symth2.r <- refit(m.symth2))@beta)
67
68            ## NB update with theta value (stored as variable, no SE) only
69            m.symth3 <- update(m.base,family=MASS::negative.binomial(theta=th1))
70            expect_equal(m.symth3@beta,(m.symth3.r <- refit(m.symth3))@beta)
71
72            ## strip NB value (off by 5e-16)
73            m.symth5 <- update(m.base,family=MASS::negative.binomial(theta=th2))
74            expect_equal(m.symth5@beta,(m.symth5.r <- refit(m.symth5))@beta)
75        })
76
77## GH #399
78test_that("na_exclude", {
79    dd1 <- dd[1:200,]
80    dd1$f1[1:5] <- NA
81    expect_error(glmer.nb(y ~ f1 + (1|g), data=dd1, na.action=na.fail),
82                 "missing values in object")
83    m1 <- glmer.nb(y ~ f1 + (1|g), data=dd1, na.action=na.omit)
84    m2 <- glmer.nb(y ~ f1 + (1|g), data=dd1, na.action=na.exclude)
85    expect_equal(fixef(m1),fixef(m1))
86    expect_equal(length(predict(m2))-length(predict(m1)),5)
87})
88
89## GH 423
90test_that("start_vals", {
91    dd1 <- dd[1:200,]
92    g1 <- glmer.nb(y ~ f1 + (1|g), data=dd1)
93    g2 <- glmer.nb(y ~ f1 + (1|g), data=dd1,
94                   initCtrl=list(theta=getME(g1,"glmer.nb.theta")))
95    expect_equal(fixef(g1),fixef(g2),tol=1e-5)
96})
97
98test_that("control arguments", {
99    dd1 <- dd[1:200,]
100    g1 <- glmer.nb(y ~ f1 + (1|g), data=dd1, initCtrl=list(theta=10))
101    expect_is(g1,"merMod")  ## dumb test - just checking for run w/o error
102    suppressWarnings(g1 <- glmer.nb(y ~ f1 + (1|g), data=dd1,
103                                    nb.control=glmerControl(optimizer="bobyqa")))
104    expect_equal(g1@optinfo$optimizer, "bobyqa")
105    suppressWarnings(g1 <- glmer.nb(y ~ f1 + (1|g), data=dd1,
106                                    nb.control=glmerControl(optCtrl=list(maxfun=2))))
107    expect_equal(g1@optinfo$feval,3)
108})
109} ## testLevel>2
110} ## testLevel>1
111