1library("testthat")
2library("lme4")
3context("specifying starting values")
4
5##' Update 'mod', copying .@call and attr(.@frame, "start")  from 'from'
6copysome <- function(mod, from) {
7    stopifnot(all.equal(class(mod), class(from)), isS4(mod))
8    mod@call <- from@call
9    attr(mod@frame, "start") <- attr(from@frame, "start")
10    mod
11}
12
13## is "Nelder_Mead" default optimizer?
14isNM <- formals(lmerControl)$optimizer == "Nelder_Mead"
15
16stMsg <- "'start' must be .* a numeric vector .* list"
17test_that("lmer", {
18    frm <- Reaction ~ Days + (Days|Subject)
19    ctrl <- lmerControl(optCtrl = list(maxfun= if(isNM) 50 else 100))
20    x <- suppressWarnings(lmer(frm, data=sleepstudy, control=ctrl, REML=FALSE))
21    x2 <- suppressWarnings(update(x,start=c(1,0,1)))
22    x3 <- suppressWarnings(update(x,start=list(theta=c(1,0,1))))
23    ff <- update(x,devFunOnly=TRUE)
24    x2@call <- x3@call <- x@call  ## hack call component
25    expect_equal(x,x2)
26    expect_equal(x,x3)
27    ## warning on deprecated list ...
28    suppressWarnings(expect_error(update(x, start = "a"), stMsg))
29    ## misspelled
30    suppressWarnings(
31        expect_error(update(x,start=list(Theta=c(1,0,1))),"incorrect components")
32    )
33    th0 <- getME(x,"theta")
34    y <- suppressWarnings(update(x,start=th0))
35    if(isNM) {
36        expect_equal(AIC(x), 1768.025, tolerance=1e-6)
37        expect_equal(AIC(y), 1763.949, tolerance=1e-6)
38    } else { ## only "bobyqa" tested:
39        expect_equal(AIC(x), 1763.939344, tolerance=1e-6)
40        expect_equal(AIC(x), AIC(y))
41    }
42    if(isNM)
43        expect_equal(suppressWarnings(optimizeLmer(ff,control=list(maxfun=50),
44                                                   start=c(1,0,1))$fval),
45                 unname(deviance(x)))
46    expect_equal(suppressWarnings(optimizeLmer(ff,control=list(maxfun=50),
47                                               start=th0)$fval),
48                 unname(deviance(y)))
49})
50test_that("glmer", {
51    ctrl <- glmerControl(optCtrl=list(maxfun=50)) # -> non-convergence warnings
52    x <- suppressWarnings(glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
53                   data = cbpp, family = binomial, control=ctrl))
54    ## theta only
55    x2 <- suppressWarnings(update(x, start= 1))
56    x3 <- suppressWarnings(update(x, start= list(theta = 1)))
57    ff <- update(x,devFunOnly=TRUE)
58    x2@call <- x3@call <- x@call  ## hack call component
59    expect_equal(x,x2)
60    expect_equal(x,x3)
61    expect_error(update(x, start="a"), stMsg)
62    expect_error(update(x, start=list(Theta=1)), "bad name\\(s\\)")
63    th0 <- getME(x,"theta")
64    y <- suppressWarnings(update(x, start=th0)) # expect_equal() fails: optinfo -> derivs -> Hessian
65
66    ## theta and beta
67    x0 <- update(x,nAGQ=0)
68    x4 <- suppressWarnings(update(x, start = list(theta=1, fixef=fixef(x0))))
69    expect_equal(x, copysome(x4, from=x))
70    x5 <- suppressWarnings(update(x, start = list(theta=1, fixef=rep(0,4))))
71    expect_equal(AIC(x5), 221.5823, tolerance=1e-6)
72    x6 <- expect_error(update(x, start = list(theta=1, fixef=rep(0,5))),
73                       "incorrect number of fixef components")
74    ## beta only
75    x7 <- suppressWarnings(update(x,start=list(fixef=fixef(x0))))
76    expect_equal(x, copysome(x7, from=x))
77    x8 <- suppressWarnings(update(x,start=list(fixef=rep(0,4))))
78    expect_equal(x5, copysome(x8, from=x5))
79})
80