[R-sig-ME] Contribution to lme4Eigen

[Joehanes, Roby (NIH/NHLBI) [F]]

Ouch. Apparently, the mail daemon ate my attachment.

So, here it is (pasted). Sorry for the mail bomb.

-----------------------
lmer <- function(formula, data, REML = TRUE, sparseX = FALSE,
control = list(), start = NULL,
verbose = 0L, subset, weights, na.action, offset,
contrasts = NULL, devFunOnly=FALSE, optimizer=c("NelderMead","bobyqa"), ...)
{
    if (sparseX) warning("sparseX = TRUE has no effect at present")
    mf <- mc <- match.call()
    ## '...' handling up front, safe-guarding against typos ("familiy") :
    if(length(l... <- list(...))) {
		if (!is.null(l...$family)) {  # call glmer if family specified
			mc[[1]] <- as.name("glmer")
			return(eval(mc, parent.frame()) )
		}
		## Check for method argument which is no longer used
		if (!is.null(method <- l...$method)) {
			msg <- paste("Argument", sQuote("method"), "is deprecated.")
			if (match.arg(method, c("Laplace", "AGQ")) == "Laplace") {
				warning(msg)
				l... <- l...[names(l...) != "method"]
			} else stop(msg)
		}
		if(length(l...))
	    warning("extra argument(s) ",
		paste(sQuote(names(l...)), collapse=", "),
		" disregarded")
    }
	
    stopifnot(length(formula <- as.formula(formula)) == 3)
    if (missing(data)) data <- environment(formula)
	# evaluate and install the model frame :
    m <- match(c("data", "subset", "weights", "na.action", "offset"),
	names(mf), 0)
    mf <- mf[c(1, m)]
    mf$drop.unused.levels <- TRUE
    mf[[1]] <- as.name("model.frame")
    fr.form <- subbars(formula) # substituted "|" by "+" -
    environment(fr.form) <- environment(formula)
    mf$formula <- fr.form
    fr <- eval(mf, parent.frame())
	# random effects and terms modules
    reTrms <- mkReTrms(findbars(formula[[3]]), fr)
    if (any(unlist(lapply(reTrms$flist, nlevels)) >= nrow(fr)))
	stop("number of levels of each grouping factor must be less than number of obs")
    ## fixed-effects model matrix X - remove random effects from formula:
    form <- formula
    form[[3]] <- if(is.null(nb <- nobars(form[[3]]))) 1 else nb
    X <- model.matrix(form, fr, contrasts)#, sparse = FALSE, row.names = FALSE) ## sparseX not yet
    p <- ncol(X)
    if ((qrX <- qr(X))$rank < p)
	stop(gettextf("rank of X = %d < ncol(X) = %d", qrX$rank, p))
    rho <- new.env(parent=parent.env(environment()))
    rho$pp <- do.call(merPredD$new, c(reTrms[c("Zt","theta","Lambdat","Lind")], n=nrow(X), list(X=X)))
    rho$resp <- mkRespMod(fr, if(REML) p else 0L)
	
    devfun <- mkdevfun(rho, 0L)
    devfun(reTrms$theta) # one evaluation to ensure all values are set
	
    if (devFunOnly) return(devfun)
	
    lower <- reTrms$lower
	# RJ's changes begin
	opt <- switch(match.arg(optimizer),
		bobyqa = {
			if(!is.numeric(control$rhobeg)) control$rhobeg <- 0.0002
			if(!is.numeric(control$rhoend)) control$rhoend <- 2e-7
			rho$control <- control
			# Delete unused options to prevent warning from showing up
			control$FtolAbs <- NULL
			control$FtolRel <- NULL
			bobyqa(rho$pp$theta, devfun, lower, control=control)
		},
		NelderMead = {
			## FIXME: this code is replicated in lmer/glmer/nlmer ...
			## it seems good to have it in R rather than C++ code but maybe it should go within Nelder_Mead() ??
			control$iprint <- switch(as.character(min(verbose,3L)), "0"=0, "1"=20,"2"=10,"3"=1)
			xst <- rep.int(0.1, length(lower))
			# RJ -- allow user control of xst, xt
			ctl <- control
			if(!is.numeric(control$xstFactor))
				xstFactor <- 0.2
			else
				xstFactor <- control$xstFactor
			if(!is.numeric(control$xtFactor))
				xtFactor <- 0.0001
			else
				xtFactor <- control$xtFactor
			ctl$xstFactor <- NULL
			ctl$xtFactor <- NULL
			Nelder_Mead(devfun, x0=rho$pp$theta, xst=xstFactor*xst, xt=xst*xtFactor, lower=lower, control=ctl)
		})

	if(optimizer=="NelderMead") {
	    if (opt$ierr < 0L) {
			if (opt$ierr > -4L)
			stop("convergence failure, code ", opt$ierr, " in NelderMead")
			else
			warning("failure to converge in", opt$control$maxfun, "evaluations")
		}
	} else if (optimizer=="bobyqa") {
		if (opt$ierr > 0L)
			warning("convergence problem, code ", opt$ierr, " in bobyqa")
	}
	# RJ's changes end
    mkMerMod(environment(devfun), opt, reTrms, fr, mc)
}## { lmer }



On Mar 1, 2012, at 2:11 PM, Joehanes, Roby (NIH/NHLBI) [F] wrote:

> Hi lme4 developers,
> 
> I hope I can contribute to something positive to the community. I noticed that even in the latest version of lme4Eigen the lmer function only use Nelder-Mead optimizer. I would like to get Bobyqa optimizer back as an option (because I really like Bobyqa optimizer). So, I patched the code a little bit. I used the (hopefully) latest revision from SVN version 1631. I also added control for xst and xt factor multiplier (a FIXME item list). I hope this change is also acceptable. If you feel I am ignorant of your code style, I apologize.
> 
> Here is the code (only the lmer function). I clearly marked my changes. If you need a diff file with the current HEAD, I'll be happy to provide you with one. The function seems to be running happily without error or warning on my machine. I hope the changes can make it to the main trunk.
> 
> Sincerely,
> Roby