### Eclipse Workspace Patch 1.0 #P lme4 Index: R/lmer.R =================================================================== --- R/lmer.R (revision 1717) +++ R/lmer.R (working copy) @@ -78,7 +78,7 @@ lmer <- function(formula, data, REML = TRUE, sparseX = FALSE, control = list(), start = NULL, verbose = 0L, subset, weights, na.action, offset, - contrasts = NULL, devFunOnly=FALSE, + contrasts = NULL, devFunOnly=FALSE, doFit = TRUE, checkLevels = FALSE, optimizer="Nelder_Mead", ...) { if (sparseX) warning("sparseX = TRUE has no effect at present") @@ -117,7 +117,7 @@ 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))) + if (checkLevels & 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 @@ -133,11 +133,12 @@ devfun <- mkdevfun(rho, 0L) devfun(reTrms$theta) # one evaluation to ensure all values are set + if (!doFit) return(list(rho=rho, devfun=devfun, reTrms=reTrms, optimizer=optimizer, control=control, fr=fr, mc=mc, nAGQ=-1, verbose=verbose)) if (devFunOnly) return(devfun) opt <- optwrap(optimizer, devfun, rho$pp$theta, lower=reTrms$lower, control=control, - rho=rho, adj=FALSE) + rho=rho, adj=FALSE, verbose=verbose) mkMerMod(environment(devfun), opt, reTrms, fr, mc) }## { lmer } @@ -262,7 +263,7 @@ glmer <- function(formula, data, family = gaussian, sparseX = FALSE, control = list(), start = NULL, verbose = 0L, nAGQ = 1L, compDev = TRUE, subset, weights, na.action, offset, - contrasts = NULL, mustart, etastart, devFunOnly = FALSE, + contrasts = NULL, mustart, etastart, devFunOnly = FALSE, doFit = TRUE, tolPwrss = 1e-7, optimizer=c("bobyqa","Nelder_Mead"), ...) { verbose <- as.integer(verbose) @@ -345,6 +346,7 @@ if (length(optimizer)==1) { optimizer <- replicate(2,optimizer) } + if (!doFit) return(list(rho=rho, devfun=devfun, reTrms=reTrms, optimizer=optimizer, control=control, fr=fr, mc=mc, nAGQ=nAGQ, family=family, verbose=verbose)) opt <- optwrap(optimizer[[1]],devfun,rho$pp$theta, rho$lower, control=control, rho=rho, adj=FALSE, verbose=verbose) @@ -405,7 +407,7 @@ ##' @export nlmer <- function(formula, data, control = list(), start = NULL, verbose = 0L, nAGQ = 1L, subset, weights, na.action, offset, - contrasts = NULL, devFunOnly = 0L, tolPwrss = 1e-10, + contrasts = NULL, devFunOnly = 0L, tolPwrss = 1e-10, doFit = TRUE, optimizer="Nelder_Mead", ...) { vals <- nlformula(mc <- match.call()) @@ -433,6 +435,7 @@ if (length(optimizer)==1) { optimizer <- replicate(2,optimizer) } + if (!doFit) return(list(rho=rho, devfun=devfun, reTrms=reTrms, optimizer=optimizer, control=control, fr=fr, mc=mc, nAGQ=nAGQ, verbose=verbose)) opt <- optwrap(optimizer[[1]],devfun, rho$pp$theta, rho$lower, control=control, rho=rho, adj=FALSE) @@ -459,6 +462,70 @@ mkMerMod(environment(devfun), opt, vals$reTrms, vals$frame, mc) }## {nlmer} +mer_finalize = function(rho, optimizer, devfun, reTrms, control, fr, mc, nAGQ, verbose) +{ + if (is(rho$resp, "lmerResp")) + { + opt <- optwrap(optimizer, + devfun, rho$pp$theta, lower=reTrms$lower, control=control, + rho=rho, adj=FALSE, verbose=verbose) + + mkMerMod(environment(devfun), opt, reTrms, fr, mc) + } else if (is(rho$resp, "glmResp")) + { + opt <- optwrap(optimizer[[1]],devfun,rho$pp$theta, rho$lower, + control=control, rho=rho, + adj=FALSE, verbose=verbose) + if (nAGQ > 0L) { + rho$lower <- c(rho$lower, rep.int(-Inf, length(rho$pp$beta0))) + rho$lp0 <- rho$pp$linPred(1) + rho$dpars <- seq_along(rho$pp$theta) + rho$baseOffset <- rho$resp$offset + 0 # forcing a copy (!) + rho$GQmat <- GHrule(nAGQ) + rho$fac <- reTrms$flist[[1]] + if (nAGQ > 1L) { + if (length(reTrms$flist) != 1L || length(reTrms$cnms[[1]]) != 1L) + stop("nAGQ > 1 is only available for models with a single, scalar random-effects term") + } + devfun <- function(pars) { + resp$updateMu(lp0) + pp$setTheta(as.double(pars[dpars])) # theta is first part of pars + resp$setOffset(baseOffset + pp$X %*% as.numeric(pars[-dpars])) + pwrssUpdate(pp, resp, tolPwrss, GQmat, compDev, fac) + } + environment(devfun) <- rho + if (devFunOnly) return(devfun) + + opt <- optwrap(optimizer[[2]],devfun,c(rho$pp$theta, rho$pp$delb), + rho$lower, control=control, rho=rho, + adj=TRUE, verbose=verbose) + } + mkMerMod(environment(devfun), opt, reTrms, fr, mc) + } else if (is(rho$resp, "nlsResp")) + { + opt <- optwrap(optimizer[[1]],devfun, rho$pp$theta, rho$lower, + control=control, rho=rho, adj=FALSE) + if (nAGQ > 0L) { + rho$lower <- c(rho$lower, rep.int(-Inf, length(rho$beta0))) + rho$u0 <- rho$pp$u0 + rho$beta0 <- rho$pp$beta0 + rho$dpars <- seq_along(rho$pp$theta) + if (nAGQ > 1L) { + if (length(vals$reTrms$flist) != 1L || length(vals$reTrms$cnms[[1]]) != 1L) + stop("nAGQ > 1 is only available for models with a single, scalar random-effects term") + rho$fac <- vals$reTrms$flist[[1]] + } + devfun <- mkdevfun(rho, nAGQ) + if (devFunOnly) return(devfun) + + opt <- optwrap(optimizer[[2]],devfun, par=c(rho$pp$theta, rho$beta0), + lower=rho$lower, control=control, rho=rho, + adj=TRUE, verbose=verbose) + } + mkMerMod(environment(devfun), opt, vals$reTrms, vals$frame, mc) + } +}## { mer_finalize } + ##' Create a deviance evaluation function from a predictor and a response module ##' ##' From an merMod object create an R function that takes a single argument, @@ -1120,15 +1187,25 @@ x0 <- c(x0, pp$beta0) lower <- c(lower, rep(-Inf,length(x0)-length(lower))) } - control <- list(...)$control - if (is.null(control)) control <- list() - control <- c(control,list(xst=0.2*xst, xt=xst*0.0001)) + l... <- list(...) + optimizer <- l...$optimizer + if (is.null(optimizer)) optimizer <- "Nelder_Mead" + if ((is(object@resp, "glmResp") | is(object@resp, "nlsResp")) & length(optimizer)==1) { + optimizer <- replicate(2,optimizer) + } + verbose <- l...$verbose + if (is.null(verbose)) verbose <- 0L + control <- l...$control + if (is.null(control)) { + control <- list(xst=0.2*xst, xt=xst*0.0001) + } ## FIXME: generic optimizer stuff ### FIXME: Probably should save the control settings and the optimizer name in the merMod object - opt <- Nelder_Mead(ff, x0, lower=lower, control=control) - mkMerMod(environment(ff), opt, list(flist=object@flist, cnms=object@cnms, Gp=object@Gp, - lower=object@lower), - object@frame, getCall(object)) +# opt <- Nelder_Mead(ff, x0, lower=lower, control=control) + opt <- optwrap(optimizer, ff, x0, lower=lower, control=control, verbose=verbose) + + mkMerMod(environment(ff), opt, list(flist=object@flist, cnms=object@cnms, Gp=object@Gp, + lower=object@lower), object@frame, getCall(object)) } ##' @S3method refitML merMod