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