## In both of these cases boundary fit (i.e. estimate of zero RE ## variance) is *incorrect*. (Nelder_Mead, restart_edge=FALSE) is the ## only case where we get stuck; either optimizer=bobyqa or ## restart_edge=TRUE (default) works if (.Platform$OS.type != "windows") { library(lme4) library(testthat) if(!dev.interactive(orNone=TRUE)) pdf("boundary_plots.pdf") ## Stephane Laurent: dat <- read.csv(system.file("testdata","dat20101314.csv", package="lme4")) fit <- lmer(y ~ (1|Operator)+(1|Part)+(1|Part:Operator), data=dat, control= lmerControl(optimizer="Nelder_Mead")) fit_b <- lmer(y ~ (1|Operator)+(1|Part)+(1|Part:Operator), data=dat, control= lmerControl(optimizer="bobyqa", restart_edge=FALSE)) fit_c <- lmer(y ~ (1|Operator)+(1|Part)+(1|Part:Operator), data=dat, control= lmerControl(optimizer="Nelder_Mead", restart_edge=FALSE, check.conv.hess="ignore")) ## final fit gives degenerate-Hessian warning ## FIXME: use fit_c with expect_warning() as a check on convergence tests ## tolerance=1e-5 seems OK in interactive use but not in R CMD check ... ?? stopifnot(all.equal(getME(fit, "theta") -> th.f, getME(fit_b,"theta"), tolerance=5e-5), all(th.f > 0)) ## Manuel Koller source(system.file("testdata", "koller-data.R", package="lme4")) ldata <- getData(13) ## old (backward compatible/buggy) fm4 <- lmer(y ~ (1|Var2), ldata, control=lmerControl(optimizer="Nelder_Mead", use.last.params=TRUE), start=list(theta=1)) fm4b <- lmer(y ~ (1|Var2), ldata, control = lmerControl(optimizer="Nelder_Mead", use.last.params=TRUE, restart_edge = FALSE, check.conv.hess="ignore", check.conv.grad="ignore"), start = list(theta=1)) ## FIXME: use as convergence test check stopifnot(getME(fm4b,"theta") == 0) fm4c <- lmer(y ~ (1|Var2), ldata, control=lmerControl(optimizer="bobyqa", use.last.params=TRUE), start=list(theta=1)) stopifnot(all.equal(getME(fm4, "theta") -> th4, getME(fm4c,"theta"), tolerance=1e-4), th4 > 0) ## new: doesn't get stuck at edge any more, but gets stuck somewhere else ... fm5 <- lmer(y ~ (1|Var2), ldata, control=lmerControl(optimizer="Nelder_Mead", check.conv.hess="ignore", check.conv.grad="ignore"), start=list(theta=1)) fm5b <- lmer(y ~ (1|Var2), ldata, control=lmerControl(optimizer="Nelder_Mead", restart_edge=FALSE, check.conv.hess="ignore", check.conv.grad="ignore"), start = list(theta = 1)) fm5c <- lmer(y ~ (1|Var2), ldata, control=lmerControl(optimizer="bobyqa"), start = list(theta = 1)) stopifnot(all.equal(unname(getME(fm5c,"theta")), 0.21067645, tolerance = 1e-7)) # 0.21067644264 [64-bit, lynne] if (require("optimx")) { ## additional stuff for diagnosing Nelder-Mead problems. fm5d <- update(fm5,control=lmerControl(optimizer="optimx", optCtrl=list(method="L-BFGS-B"))) fm5e <- update(fm5, control=lmerControl(optimizer="nloptwrap")) mList <- setNames(list(fm4,fm4b,fm4c,fm5,fm5b,fm5c,fm5d,fm5e), c("NM/uselast","NM/uselast/norestart","bobyqa/uselast", "NM","NM/norestart","bobyqa","LBFGSB","nloptr/bobyqa")) pp <- profile(fm5c,which=1) dd <- as.data.frame(pp) par(las=1,bty="l") v <- sapply(mList, function(x) sqrt(VarCorr(x)[[1]])) plot(.zeta^2~.sig01, data=dd, type="b") abline(v=v) res <- cbind(VCorr = sapply(mList, function(x) sqrt(VarCorr(x)[[1]])), theta = sapply(mList, getME,"theta"), loglik = sapply(mList, logLik)) res print(sessionInfo(), locale=FALSE) } ###################### library(lattice) ## testing boundary and near-boundary cases tmpf <- function(i,...) { set.seed(i) d <- data.frame(x=rnorm(60),f=factor(rep(1:6,each=10))) d$y <- simulate(~x+(1|f),family=gaussian,newdata=d, newparams=list(theta=0.01,beta=c(1,1),sigma=5))[[1]] lmer(y~x+(1|f),data=d,...) } sumf <- function(m) { unlist(VarCorr(m))[1] } if (FALSE) { ## figuring out which seeds will give boundary and ## near-boundary solutions mList <- lapply(1:201,tmpf) # [FIXME tons of messages "theta parameters vector not named"] ss <- sapply(mList,sumf)+1e-50 par(las=1,bty="l") hist(log(ss),col="gray",breaks=50) ## values lying on boundary which(log(ss)<(-40)) ## 5, 7-13, 15, 21, ... ## values close to boundary (if check.edge not set) which(log(ss)>(-40) & log(ss) <(-20)) ## 16, 44, 80, 86, 116, ... } ## diagnostic plot tmpplot <- function(i, FUN=tmpf) { dd <- FUN(i, devFunOnly=TRUE) x <- 10^seq(-10,-6.5,length=201) dvec <- sapply(x,dd) op <- par(las=1,bty="l"); on.exit(par(op)) plot(x,dvec-min(dvec)+1e-16, log="xy", type="b") r <- FUN(i) abline(v = getME(r,"theta"), col=2) invisible(r) } ## Case #1: boundary estimate with or without boundary.tol m5 <- tmpf(5) m5B <- tmpf(5,control=lmerControl(boundary.tol=0)) stopifnot(getME(m5, "theta")==0, getME(m5B,"theta")==0) p5 <- profile(m5) ## bobyqa warnings but results look reasonable xyplot(p5) ## reveals slight glitch (bottom row of plots doesn't look right) expect_warning(splom(p5),"unreliable for singular fits") p5B <- profile(m5, signames=FALSE) # -> bobyqa convergence warning (code 3) expect_warning(splom(p5B), "unreliable for singular fits") if(lme4:::testLevel() >= 2) { ## avoid failure to warn ## Case #2: near-boundary estimate, but boundary.tol can't fix it m16 <- tmpplot(16) ## sometimes[2014-11-11] fails (??) : p16 <- profile(m16) ## warning message*s* (non-monotonic profile and more) plotOb <- xyplot(p16) ## NB: It's the print()ing of 'plotOb' which warns ==> need to do this explicitly: expect_warning(print(plotOb), ## warns about linear interpolation in profile for variable 1 "using linear interpolation") d16 <- as.data.frame(p16) xyplot(.zeta ~ .focal|.par, data=d16, type=c("p","l"), scales = list(x=list(relation="free"))) try(splom(p16)) ## breaks when calling predict(.) } ## bottom line: ## * xyplot.thpr could still be improved ## * most of the near-boundary cases are noisy and can't easily be ## fixed tmpf2 <- function(i,...) { set.seed(i) d <- data.frame(x=rnorm(60),f=factor(rep(1:6,each=10)), w=rep(10,60)) d$y <- simulate(~x+(1|f),family=binomial, weights=d$w,newdata=d, newparams=list(theta=0.01,beta=c(1,1)))[[1]] glmer(y~x+(1|f),data=d,family=binomial,weights=w,...) } if (FALSE) { ## figuring out which seeds will give boundary and ## near-boundary solutions mList <- lapply(1:201,tmpf2) ss <- sapply(mList,sumf)+1e-50 par(las=1,bty="l") hist(log(ss),col="gray",breaks=50) ## values lying on boundary head(which(log(ss)<(-50))) ## 1-5, 7 ... ## values close to boundary (if check.edge not set) which(log(ss)>(-50) & log(ss) <(-20)) ## 44, 46, 52, ... } ## m1 <- tmpf2(1) ## FIXME: doesn't work if we generate m1 via tmpf2(1) -- ## some environment lookup problem ... set.seed(1) d <- data.frame(x=rnorm(60),f=factor(rep(1:6,each=10)), w=rep(10,60)) d$y <- simulate(~x+(1|f),family=binomial, weights=d$w,newdata=d, newparams=list(theta=0.01,beta=c(1,1)))[[1]] m1 <- glmer(y~x+(1|f),data=d,family=binomial,weights=w) p1 <- profile(m1) xyplot(p1) expect_warning(splom(p1),"splom is unreliable") } ## skip on windows (for speed)