Nothing
## 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)
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.