Nothing
library("testthat")
library("lme4")
context("specifying starting values")
##' Update 'mod', copying .@call and attr(.@frame, "start") from 'from'
copysome <- function(mod, from) {
stopifnot(all.equal(class(mod), class(from)), isS4(mod))
mod@call <- from@call
attr(mod@frame, "start") <- attr(from@frame, "start")
mod
}
## is "Nelder_Mead" default optimizer?
isNM <- formals(lmerControl)$optimizer == "Nelder_Mead"
stMsg <- "'start' must be .* a numeric vector .* list"
test_that("lmer", {
frm <- Reaction ~ Days + (Days|Subject)
ctrl <- lmerControl(optCtrl = list(maxfun= if(isNM) 50 else 100))
x <- suppressWarnings(lmer(frm, data=sleepstudy, control=ctrl, REML=FALSE))
x2 <- suppressWarnings(update(x,start=c(1,0,1)))
x3 <- suppressWarnings(update(x,start=list(theta=c(1,0,1))))
ff <- update(x,devFunOnly=TRUE)
x2@call <- x3@call <- x@call ## hack call component
expect_equal(x,x2)
expect_equal(x,x3)
## warning on deprecated list ...
suppressWarnings(expect_error(update(x, start = "a"), stMsg))
## misspelled
suppressWarnings(
expect_error(update(x,start=list(Theta=c(1,0,1))),"incorrect components")
)
th0 <- getME(x,"theta")
y <- suppressWarnings(update(x,start=th0))
if(isNM) {
expect_equal(AIC(x), 1768.025, tolerance=1e-6)
expect_equal(AIC(y), 1763.949, tolerance=1e-6)
} else { ## only "bobyqa" tested:
expect_equal(AIC(x), 1763.939344, tolerance=1e-6)
expect_equal(AIC(x), AIC(y))
}
if(isNM)
expect_equal(suppressWarnings(optimizeLmer(ff,control=list(maxfun=50),
start=c(1,0,1))$fval),
unname(deviance(x)))
expect_equal(suppressWarnings(optimizeLmer(ff,control=list(maxfun=50),
start=th0)$fval),
unname(deviance(y)))
})
test_that("glmer", {
ctrl <- glmerControl(optCtrl=list(maxfun=50)) # -> non-convergence warnings
x <- suppressWarnings(glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
data = cbpp, family = binomial, control=ctrl))
## theta only
x2 <- suppressWarnings(update(x, start= 1))
x3 <- suppressWarnings(update(x, start= list(theta = 1)))
ff <- update(x,devFunOnly=TRUE)
x2@call <- x3@call <- x@call ## hack call component
expect_equal(x,x2)
expect_equal(x,x3)
expect_error(update(x, start="a"), stMsg)
expect_error(update(x, start=list(Theta=1)), "bad name\\(s\\)")
th0 <- getME(x,"theta")
y <- suppressWarnings(update(x, start=th0)) # expect_equal() fails: optinfo -> derivs -> Hessian
## theta and beta
x0 <- update(x,nAGQ=0)
x4 <- suppressWarnings(update(x, start = list(theta=1, fixef=fixef(x0))))
expect_equal(x, copysome(x4, from=x))
x5 <- suppressWarnings(update(x, start = list(theta=1, fixef=rep(0,4))))
expect_equal(AIC(x5), 221.5823, tolerance=1e-6)
x6 <- expect_error(suppressWarnings(update(x, start = list(theta=1, fixef=rep(0,5))),
"incorrect number of fixef components"))
## beta only
x7 <- suppressWarnings(update(x,start=list(fixef=fixef(x0))))
expect_equal(x, copysome(x7, from=x))
x8 <- suppressWarnings(update(x,start=list(fixef=rep(0,4))))
expect_equal(x5, copysome(x8, from=x5))
})
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.