lme4: Linear Mixed-Effects Models using 'Eigen' and S4

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