[Lme4-commits] r1644 - in pkg/lme4Eigen: . R man tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Mar 11 03:18:32 CET 2012
Author: bbolker
Date: 2012-03-11 03:18:31 +0100 (Sun, 11 Mar 2012)
New Revision: 1644
Added:
pkg/lme4Eigen/R/plot.R
Modified:
pkg/lme4Eigen/DESCRIPTION
pkg/lme4Eigen/NAMESPACE
pkg/lme4Eigen/R/bootMer.R
pkg/lme4Eigen/R/lmer.R
pkg/lme4Eigen/R/optimizer.R
pkg/lme4Eigen/R/utilities.R
pkg/lme4Eigen/do-roxy
pkg/lme4Eigen/man/Nelder_Mead.Rd
pkg/lme4Eigen/man/getME.Rd
pkg/lme4Eigen/man/glmer.Rd
pkg/lme4Eigen/man/lmer.Rd
pkg/lme4Eigen/man/mkdevfun.Rd
pkg/lme4Eigen/man/nlmer.Rd
pkg/lme4Eigen/tests/drop.R
pkg/lme4Eigen/tests/glmer-1.R
pkg/lme4Eigen/tests/nbinom.R
Log:
restructuring of optimization
added getOptfun: convert from string to function (if necessary),
check function for basic properties (correct
formal arguments)
added optwrap: wrapper function for optimization settings
Nelder_Mead:
* move verbose -> iprint logic inside Nelder_Mead function
* move xst, xt settings inside function; add xst, xt to control()
* changed argument and return value names for conformity with API:
fn, par, lower, upper, control -> fval, par, convergence, message, control
add optimizer argument to lmer, glmer, nlmer: of length 2 for glmer, nlmer
added xyplot of data to glmer examples
modified examples(ranef) to allow ranef(.,postVar=TRUE) for single-RE case
added getME(.,"lower")
added p-value column to GLMM summary
added as.data.frame.boot (not exposed)
mkMerMod: test for conv*, fval* rather than conv, fval
add plot methods
add BMB settings to do-roxy
tweaks to some tests
add optimx to 'suggests'
add FIXME in mkdevfun.Rd
Modified: pkg/lme4Eigen/DESCRIPTION
===================================================================
--- pkg/lme4Eigen/DESCRIPTION 2012-03-10 20:25:38 UTC (rev 1643)
+++ pkg/lme4Eigen/DESCRIPTION 2012-03-11 02:18:31 UTC (rev 1644)
@@ -1,5 +1,5 @@
Package: lme4Eigen
-Version: 0.9996875-11
+Version: 0.9996875-13
Date: $Date$
Revision: $Revision$
Title: Linear mixed-effects models using Eigen and S4
@@ -34,7 +34,8 @@
mlmRev,
plyr,
reshape,
- testthat
+ testthat,
+ optimx
LazyData: yes
License: GPL (>=2)
URL: http://lme4.r-forge.r-project.org/
@@ -54,3 +55,5 @@
'bootMer.R'
'nbinom.R'
'mcmcsamp.R'
+ 'plot.R'
+BuildVignettes: no
Modified: pkg/lme4Eigen/NAMESPACE
===================================================================
--- pkg/lme4Eigen/NAMESPACE 2012-03-10 20:25:38 UTC (rev 1643)
+++ pkg/lme4Eigen/NAMESPACE 2012-03-11 02:18:31 UTC (rev 1644)
@@ -118,6 +118,7 @@
S3method(plot,coef.mer)
S3method(plot,lmList.confint)
S3method(plot,ranef.mer)
+S3method(plot,merMod)
S3method(predict,merMod)
S3method(print,merMod)
S3method(print,summary.mer)
Modified: pkg/lme4Eigen/R/bootMer.R
===================================================================
--- pkg/lme4Eigen/R/bootMer.R 2012-03-10 20:25:38 UTC (rev 1643)
+++ pkg/lme4Eigen/R/bootMer.R 2012-03-11 02:18:31 UTC (rev 1644)
@@ -63,7 +63,7 @@
##' fm01ML <- lmer(Yield ~ 1|Batch, Dyestuff, REML = FALSE)
##' ## see ?"profile-methods"
##' mySumm <- function(.) { s <- sigma(.)
-##' c(beta =getME(., "beta"), sigma = s, sig01 = s * getME(., "theta")) }
+##' c(beta =getME(., "beta"), sigma = s, sig01 = unname(s * getME(., "theta"))) }
##' (t0 <- mySumm(fm01ML)) # just three parameters
##'
##' ## 3.8s (on a 5600 MIPS 64bit fast(year 2009) desktop "AMD Phenom(tm) II X4 925"):
@@ -146,3 +146,8 @@
ran.gen = "simulate(<lmerMod>, 1, *)", mle = mle),
class = "boot")
} ## {bootMer}
+
+as.data.frame.boot <- function(x,...) {
+ as.data.frame(x$t)
+}
+
Modified: pkg/lme4Eigen/R/lmer.R
===================================================================
--- pkg/lme4Eigen/R/lmer.R 2012-03-10 20:25:38 UTC (rev 1643)
+++ pkg/lme4Eigen/R/lmer.R 2012-03-11 02:18:31 UTC (rev 1644)
@@ -60,6 +60,7 @@
##' @param contrasts an optional list. See the \code{contrasts.arg} of
##' \code{model.matrix.default}.
##' @param devFunOnly logical - return only the deviance evaluation function.
+##' @param optimizer character - name of optimizing function
##' @param \dots other potential arguments. A \code{method} argument was used
##' in earlier versions of the package. Its functionality has been replaced by
##' the \code{REML} argument.
@@ -77,7 +78,8 @@
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,
+ optimizer="Nelder_Mead", ...)
{
if (sparseX) warning("sparseX = TRUE has no effect at present")
mf <- mc <- match.call()
@@ -133,26 +135,13 @@
if (devFunOnly) return(devfun)
- ## 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() ??
+ opt <- optwrap(optimizer,
+ devfun, rho$pp$theta, lower=reTrms$lower, control=control, rho, adj=FALSE)
- control$iprint <- switch(as.character(min(verbose,3L)),
- "0"=0, "1"=20,"2"=10,"3"=1)
-
- lower <- reTrms$lower
- ## FIXME: allow user control of xst, xt ?
- xst <- rep.int(0.1, length(lower))
- opt <- Nelder_Mead(devfun, x0=rho$pp$theta, xst=0.2*xst, xt=xst*0.0001,
- lower=lower, control=control)
- 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")
- }
mkMerMod(environment(devfun), opt, reTrms, fr, mc)
}## { lmer }
+
##' Fit a generalized linear mixed model (GLMM)
##'
##' Fit a generalized linear mixed model, which incorporates both fixed-effects
@@ -203,11 +192,35 @@
##' \code{theta}, these are used as the starting values for those slots.
##' A numeric \code{start} argument of the appropriate length is used as the
##' starting value of \code{theta}.
-##' @param optimizer which optimizer to use in the second phase of optimization.
-##' The default is \code{\link{NelderMead}} and the alternative is
-##' \code{\link{bobyqa}}. For difficult model fits we have found
-##' \code{\link{NelderMead}} to be more reliable but occasionally slower than
-##' \code{\link{bobyqa}}.
+##' @param optimizer which optimizer(s) to use for each phase of optimization.
+##' A character vector or list of functions.
+##' If \code{length(optimizer)==2}, the first element will be used
+##' for the preliminary (random effects parameters only) optimization, while
+##' the second will be used for the final (random effects plus
+##' fixed effect parameters) phase. The built-in optimizers are
+##' \code{\link{Nelder_Mead}} and \code{\link{bobyqa}} (from
+##' the \code{minqa} package; the default
+##' is to use \code{\link{bobyqa}} for the first and
+##' \code{\link{Nelder_Mead}} for the final phase.
+##' (FIXME: simplify if possible!). For difficult model fits we have found
+##' \code{\link{Nelder_Mead}} to be more reliable but occasionally slower than
+##' \code{\link{bobyqa}}. Any minimizing function that allows
+##' box constraints can be used
+##' provided that it (1) takes input parameters \code{fn} (function
+##' to be optimized), \code{par} (starting parameter values),
+##' \code{lower} (lower bounds) and \code{control} (control parameters,
+##' passed through from the \code{control} argument) and (2)
+##' returns a list with (at least) elements \code{par} (best-fit
+##' parameters), \code{fval} (best-fit function value), \code{conv}
+##' (convergence code) and (optionally) \code{message} (informational
+##' message, or explanation of convergence failure).
+##' Special provisions are made for \code{\link{bobyqa}},
+##' \code{\link{Nelder_Mead}}, and optimizers wrapped in
+##' the \code{optimx} package; to use \code{optimx} optimizers
+##' (including \code{L-BFGS-B} from base \code{optim} and
+##' \code{nlminb}), pass the \code{method} argument to \code{optim}
+##' in the \code{control} argument.
+##'
##' @param mustart optional starting values on the scale of the conditional mean,
##' as in \code{\link[stats]{glm}}; see there for details.
##' @param etastart optional starting values on the scale of the unbounded
@@ -225,6 +238,8 @@
##' @keywords models
##' @examples
##' ## generalized linear mixed model
+##' library(lattice)
+##' xyplot(incidence/size ~ period, group=herd, type="a", data=cbpp)
##' (gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
##' data = cbpp, family = binomial))
##' ## using nAGQ=0L only gets close to the optimum
@@ -248,7 +263,7 @@
control = list(), start = NULL, verbose = 0L, nAGQ = 1L,
compDev = TRUE, subset, weights, na.action, offset,
contrasts = NULL, mustart, etastart, devFunOnly = FALSE,
- tolPwrss = 1e-10, optimizer=c("NelderMead","bobyqa"), ...)
+ tolPwrss = 1e-10, optimizer=c("bobyqa","Nelder_Mead"), ...)
{
verbose <- as.integer(verbose)
mf <- mc <- match.call()
@@ -329,11 +344,12 @@
devfun <- mkdevfun(rho, 0L) # deviance as a function of theta only
if (devFunOnly && !nAGQ) return(devfun)
-
- ## FIXME: why is bobyqa always used for preliminary fit? document??
- control$iprint <- min(verbose, 3L)
-
- opt <- bobyqa(rho$pp$theta, devfun, rho$lower, control=control)
+ if (length(optimizer)==1) {
+ optimizer <- replicate(2,optimizer)
+ }
+ opt <- optwrap(optimizer[[1]],devfun,rho$pp$theta, rho$lower, control, rho,
+ adj=FALSE, verbose=verbose)
+
rho$nAGQ <- nAGQ
if (nAGQ > 0L) {
rho$lower <- c(rho$lower, rep.int(-Inf, length(rho$beta0)))
@@ -347,23 +363,11 @@
}
devfun <- mkdevfun(rho, nAGQ)
if (devFunOnly) return(devfun)
- 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
- bobyqa(c(rho$pp$theta, rho$beta0), devfun, rho$lower, control=control)
- },
- NelderMead = {
- control$iprint <- switch(as.character(min(verbose,3L)),
- "0"=0,"1"=20,"2"=10,"3"=1)
- xst <- c(rep.int(0.1, length(rho$dpars)),
- sqrt(diag(environment(devfun)$pp$unsc())))
- Nelder_Mead(devfun, x0=with(environment(devfun), c(pp$theta, pp$beta0)),
- xst=0.2*xst, xt= xst*0.0001, lower=rho$lower, control=control)
- })
- }
+ opt <- optwrap(optimizer[[2]],devfun,c(rho$pp$theta, rho$beta0), rho$lower, control, rho,
+ adj=TRUE, verbose=verbose)
+ }
+
mkMerMod(environment(devfun), opt, reTrms, fr, mc)
}## {glmer}
@@ -398,7 +402,7 @@
nlmer <- function(formula, data, control = list(), start = NULL, verbose = 0L,
nAGQ = 1L, subset, weights, na.action, offset,
contrasts = NULL, devFunOnly = 0L, tolPwrss = 1e-10,
- optimizer=c("NelderMead","bobyqa"), ...)
+ optimizer="Nelder_Mead", ...)
{
vals <- nlformula(mc <- match.call())
if ((qrX <- qr(X <- vals$X))$rank < (p <- ncol(X)))
@@ -422,19 +426,12 @@
rho$beta0 <- rho$pp$beta0
rho$tolPwrss <- tolPwrss # Resetting this is intentional. The initial optimization is coarse.
- control$iprint <- switch(as.character(min(verbose,3L)),
- "0"=0,"1"=20,"2"=10,"3"=1)
- lower <- rho$lower
- xst <- rep.int(0.1, length(lower))
+ if (length(optimizer)==1) {
+ optimizer <- replicate(2,optimizer)
+ }
- opt <- Nelder_Mead(devfun, x0=rho$pp$theta, xst=0.2*xst, xt=xst*0.0001,
- lower=lower, control=control)
- if (opt$ierr < 0L) {
- if (opt$ierr > -4L)
- stop("convergence failure, code ", opt$ierr, " in NelderMead")
- else
- warning("failure to converge in ", cc$maxfun, " evaluations")
- }
+ opt <- optwrap(optimizer[[1]],devfun, rho$pp$theta, rho$lower, control, rho, adj=FALSE)
+
if (nAGQ > 0L) {
rho$lower <- c(rho$lower, rep.int(-Inf, length(rho$beta0)))
rho$u0 <- rho$pp$u0
@@ -447,20 +444,12 @@
}
devfun <- mkdevfun(rho, nAGQ)
if (devFunOnly) return(devfun)
- 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
- bobyqa(c(rho$pp$theta, rho$beta0), devfun, rho$lower, control=control)
- },
- NelderMead = {
- xst <- c(rep.int(0.1, length(rho$dpars)),
- sqrt(diag(environment(devfun)$pp$unsc())))
- Nelder_Mead(devfun, x0=with(environment(devfun), c(pp$theta, pp$beta0)),
- xst=0.2*xst, xt= xst*0.0001, lower=rho$lower, control=control)
- })
- }
+
+ opt <- optwrap(optimizer[[2]],devfun, par=c(rho$pp$theta, rho$beta0), lower=rho$lower, control=control, rho,
+ adj=TRUE, verbose=verbose)
+
+
+ }
mkMerMod(environment(devfun), opt, vals$reTrms, vals$frame, mc)
}## {nlmer}
@@ -961,9 +950,9 @@
##' fm2 <- lmer(Reaction ~ Days + (1|Subject) + (0+Days|Subject), sleepstudy)
##' fm3 <- lmer(diameter ~ (1|plate) + (1|sample), Penicillin)
##' ranef(fm1)
-##' if(FALSE) { ##-- postVar=TRUE is not yet implemented -- FIXME
-##' str(rr1 <- ranef(fm1, postVar = TRUE))
+##' ##' str(rr1 <- ranef(fm1, postVar = TRUE))
##' dotplot(rr1,scales = list(x = list(relation = 'free')))[["Subject"]]
+##' if(FALSE) { ##-- postVar=TRUE is not yet implemented for multiple terms -- FIXME
##' str(ranef(fm2, postVar = TRUE))
##' }
##' op <- options(digits = 4)
@@ -1078,8 +1067,10 @@
}
control <- list(...)$control
if (is.null(control)) control <- list()
+ control <- c(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, xst=0.2*xst, xt=xst*0.0001, lower=lower, control=control)
+ opt <- Nelder_Mead(ff, x0, lower=lower, control=control)
mkMerMod(environment(ff), opt, list(flist=object at flist, cnms=object at cnms, Gp=object at Gp,
lower=object at lower), object at frame, getCall(object))
}
@@ -1449,6 +1440,7 @@
##' \item{devcomp}{a list consisting of a named numeric vector, \dQuote{cmp}, and
##' a named integer vector, \dQuote{dims}, describing the fitted model}
##' \item{offset}{model offset}
+##' \item{lower}{lower bounds on model parameters} ## FIXME: theta only?
##' }
##' @return Unspecified, as very much depending on the \code{\link{name}}.
##' @seealso \code{\link{getCall}()},
@@ -1479,7 +1471,7 @@
"RX", "RZX",
"beta", "theta",
"REML", "n_rtrms", "is_REML", "devcomp",
- "offset"))
+ "offset", "lower"))
{
if(missing(name)) stop("'name' must not be missing")
stopifnot(length(name <- as.character(name)) == 1,
@@ -1529,6 +1521,7 @@
"devcomp" = dc,
"offset" = rsp$offset,
+ "lower" = object at lower, ## FIXME: should this include -Inf values for beta?
"..foo.." =# placeholder!
stop(gettextf("'%s' is not implemented yet",
sprintf("getME(*, \"%s\")", name))),
@@ -1721,6 +1714,10 @@
if (nrow(coefs) > 0) {
coefs <- cbind(coefs, coefs[,1]/coefs[,2], deparse.level=0)
colnames(coefs)[3] <- paste(if(useSc) "t" else "z", "value")
+ if (isGLMM(object)) {
+ coefs <- cbind(coefs,2*pnorm(abs(coefs[,3]),lower.tail=FALSE))
+ colnames(coefs)[4] <- c("Pr(>|z|)")
+ }
}
mName <- paste(switch(1L + dd["GLMM"] * 2L + dd["NLMM"],
"Linear", "Nonlinear",
@@ -1900,3 +1897,71 @@
## **not** is(object at resp,"lmerResp") ?
}
+getOptfun <- function(optimizer) {
+ if (is.character(optimizer)) {
+ optfun <- try(get(optimizer),silent=TRUE)
+ } else optfun <- optimizer
+ if (inherits(optfun,"try-error")) stop("couldn't find optimizer function ",optimizer)
+ if (!is(optfun,"function")) stop("non-function specified as optimizer")
+ if (any(is.na(match(c("fn","par","lower","control"),names(formals(optfun))))))
+ stop("optimizer function must use (at least) formal parameters 'fn', 'par', 'lower', 'control'")
+ optfun
+}
+
+optwrap <- function(optimizer,fn,par,lower,control,rho,adj, verbose=0L) {
+ optfun <- getOptfun(optimizer)
+ ## special-purpose control parameter tweaks: only for second round in nlmer, glmer
+ if (adj && is.character(optimizer))
+ switch(optimizer,
+ bobyqa = {
+ if(!is.numeric(control$rhobeg)) control$rhobeg <- 0.0002
+ if(!is.numeric(control$rhoend)) control$rhoend <- 2e-7
+ rho$control <- control
+ },
+ Nelder_Mead = {
+ if (is.null(control$xst))
+ xst <- c(rep.int(0.1, length(rho$dpars)),
+ sqrt(diag(environment(fn)$pp$unsc())))
+ control$xst <- 0.2*xst
+ control$verbose <- verbose
+ if (is.null(control$xt)) control$xt <- control$xst*5e-4
+ rho$control <- control
+ })
+ arglist <- list(fn=fn, par=par, lower=lower, control=control)
+ ## optimx: must pass method in control (?) because 'method' was previously
+ ## used in lme4 to specify REML vs ML
+ ## FIXME: test -- does deparse(substitute(...)) clause work?
+ if (optimizer=="optimx" || deparse(substitute(optimizer))=="optimx") {
+ if (is.null(method <- control$method)) stop("must specify 'method' explicitly for optimx")
+ arglist$control$method <- NULL
+ arglist <- c(arglist,list(method=method))
+ }
+ ## FIXME: test! effects of multiple warnings??
+ ## may not need to catch warnings after all??
+ curWarning <- NULL
+ opt <- withCallingHandlers(do.call(optfun,arglist),
+ warning = function(w) {
+ ## browser()
+ cat("caught warning:",w$message,"\n")
+ assign("curWarning",w$message,parent.frame())
+ invokeRestart("muffleWarning")
+ })
+ ## if (!is.null(curWarning)) browser()
+ ## FIXME: set code to warn on convergence !=0
+ ## post-fit tweaking
+ if (optimizer=="bobyqa") {
+ opt$convergence <- opt$ierr
+ }
+ if (optimizer=="optimx") {
+ optr <- lapply(opt,"[[",1)[c("par","fvalues","conv")]
+ optr$message <- attr(opt,"details")[[1]]$message
+ opt <- optr
+ }
+ if (opt$conv!=0) {
+ wmsg <- paste("convergence code",opt$conv,"from",optimizer)
+ if (!is.null(opt$msg)) wmsg <- paste(wmsg,": ",opt$msg,sep="")
+ warning(wmsg)
+ }
+ opt
+}
+
Modified: pkg/lme4Eigen/R/optimizer.R
===================================================================
--- pkg/lme4Eigen/R/optimizer.R 2012-03-10 20:25:38 UTC (rev 1643)
+++ pkg/lme4Eigen/R/optimizer.R 2012-03-11 02:18:31 UTC (rev 1644)
@@ -1,23 +1,27 @@
##' Nelder-Mead optimization of parameters that may be subject to box constraints
##'
##' @title Nelder-Mead optimization
-##' @param ff a function of one numeric vector argument returning a numeric scalar
+##' @param fn a function of one numeric vector argument returning a numeric scalar
+##' @param par numeric vector of starting values for the parameters.
##' @param lower numeric vector of lower bounds - elements may be \code{-Inf}.
##' @param upper numeric vector of upper bounds - elements may be \code{Inf}.
-##' @param xst numeric vector of initial step sizes to establish the simplex -
-##' all elements must be non-zero.
-##' @param x0 numeric vector of starting values for the parameters.
-##' @param xt numeric vector of tolerances on the parameters.
##' @param control a named list of control settings. Possible settings are
##' \describe{
##' \item{iprint}{numeric scalar - frequency of printing evaluation information.
##' Defaults to 0 indicating no printing.}
-##' \item{maxfun}{numeric scalar - maximum number of function evaluations allowed.}
-##' \item{FtolAbs}{numeric scalar - absolute tolerance on change in function values}
-##' \item{FtolRel}{numeric scalar - relative tolerance on change in function values}
-##' \item{XtolRel}{numeric scalar - relative tolerance on change in parameter values}
-##' \item{MinfMax}{numeric scalar - maximum value of the minimum}
+##' \item{maxfun}{numeric scalar - maximum number of function evaluations allowed (default:10000).}
+##' \item{FtolAbs}{numeric scalar - absolute tolerance on change in function values (default: 1e-5)}
+##' \item{FtolRel}{numeric scalar - relative tolerance on change in function values (default:1e-15)}
+##' \item{XtolRel}{numeric scalar - relative tolerance on change in parameter values (default: 1e-7)}
+##' \item{MinfMax}{numeric scalar - maximum value of the minimum (default: .Machine$double.xmin)}
+##' \item{xst}{numeric vector of initial step sizes to establish the simplex -
+##' all elements must be non-zero (default: rep(0.02,length(par)))}
+##' \item{xt}{numeric vector of tolerances on the parameters (default: xst*5e-4)}
+##' \item{verbose}{numeric value: 0=no printing, 1=print every 20 evaluations,
+##' 2=print every 10 evalutions, 3=print every evaluation. Sets \sQuote{iprint},
+##' if specified, but does not override it.}
##' }
+##'
##' @return a list with 4 components
##' \item{fval}{numeric scalar - the minimum function value achieved}
##' \item{par}{numeric vector - the value of \code{x} providing the minimum}
@@ -32,16 +36,30 @@
##' \item{-1}{\code{nm_x0notfeasible}: initial x is not feasible (?)}
##' }
##' @export
-Nelder_Mead <- function(ff, x0, xst, xt, lower=rep.int(-Inf, n),
+Nelder_Mead <- function(fn, par, lower=rep.int(-Inf, n),
upper=rep.int(Inf, n), control=list()) {
- stopifnot(is.function(ff),
- length(formals(ff)) == 1L,
- (n <- length(x0 <- as.numeric(x0))) == length(lower <- as.numeric(lower)),
+ n <- length(par)
+ if (is.null(xst <- control[["xst"]])) xst <- rep.int(0.02,n)
+ if (is.null(xt <- control[["xt"]])) xt <- xst*5e-4
+
+ control[["xst"]] <- control[["xt"]] <- NULL
+
+ ## mapping between simpler 'verbose' setting (0=no printing, 1=20, 2=10, 3=1)
+ ## and internal 'iprint' control (frequency of printing)
+ if (is.null(verbose <- control[["verbose"]])) verbose <- 0
+ control[["verbose"]] <- NULL
+ if (is.null(control[["iprint"]])) {
+ control[["iprint"]] <- switch(as.character(min(as.numeric(verbose),3L)),
+ "0"=0, "1"=20,"2"=10,"3"=1)
+ }
+ stopifnot(is.function(fn),
+ length(formals(fn)) == 1L,
+ (n <- length(par <- as.numeric(par))) == length(lower <- as.numeric(lower)),
length(upper <- as.numeric(upper)) == n,
length(xst <- as.numeric(xst)) == n,
all(xst != 0),
length(xt <- as.numeric(xt)) == n)
- nM <- NelderMead$new(lower=lower, upper=upper, x0=x0, xst=xst, xt=xt)
+ nM <- NelderMead$new(lower=lower, upper=upper, x0=par, xst=xst, xt=xt)
cc <- do.call(function(iprint=0L, maxfun=10000L, FtolAbs=1e-5,
FtolRel=1e-15, XtolRel=1e-7,
MinfMax=.Machine$double.xmin, ...) {
@@ -54,6 +72,22 @@
nM$setIprint(cc$iprint)
nM$setMaxeval(cc$maxfun)
nM$setMinfMax(cc$MinfMax)
- while ((nMres <- nM$newf(ff(nM$xeval()))) == 0L) {}
- list(fval=nM$value(), par=nM$xpos(), ierr=nMres, control=cc)
+ while ((nMres <- nM$newf(fn(nM$xeval()))) == 0L) {}
+
+ cmsg <- "reached max evaluations"
+ if (nMres==-4) {
+ ## map max evals from error to warning
+ cmsg <- warning(sprintf("failure to converge in %d evaluations",cc$maxfun))
+ nMres <- 4
+ }
+
+ msgvec <- c("nm_forced","cannot generate a feasible simplex","initial x is not feasible",
+ "active","minf_max","fcvg","xcvg", ## FIXME: names (see NelderMead_newf in external.cpp)
+ cmsg)
+
+ if (nMres<0) stop(msgvec[nMres+4])
+
+ cc <- c(cc,xst=xst,xt=xt)
+ list(fval=nM$value(), par=nM$xpos(), convergence=pmin(nMres,0), message=msgvec[nMres+4],
+ control=cc)
}
Added: pkg/lme4Eigen/R/plot.R
===================================================================
--- pkg/lme4Eigen/R/plot.R (rev 0)
+++ pkg/lme4Eigen/R/plot.R 2012-03-11 02:18:31 UTC (rev 1644)
@@ -0,0 +1,352 @@
+## copied/modified from nlme
+splitFormula <-
+ ## split, on the nm call, the rhs of a formula into a list of subformulas
+ function(form, sep = "/")
+{
+ if (inherits(form, "formula") ||
+ mode(form) == "call" && form[[1]] == as.name("~"))
+ return(splitFormula(form[[length(form)]], sep = sep))
+ if (mode(form) == "call" && form[[1]] == as.name(sep))
+ return(do.call("c", lapply(as.list(form[-1]), splitFormula, sep = sep)))
+ if (mode(form) == "(") return(splitFormula(form[[2]], sep = sep))
+ if (length(form) < 1) return(NULL)
+ list(asOneSidedFormula(form))
+}
+
+allVarsRec <-
+ ## Recursive version of all.vars
+ function(object)
+{
+ if (is.list(object)) {
+ unlist(lapply(object, allVarsRec))
+ } else {
+ all.vars(object)
+ }
+}
+
+## crippled version of getData.gnls from nlme
+getData <- function(object)
+{
+ mCall <- object at call
+ data <- eval(mCall$data)
+ if (is.null(data)) return(data)
+ ## FIXME: deal with NAs, subset appropriately
+ ## naPat <- eval(mCall$naPattern)
+ ## if (!is.null(naPat)) {
+ ## data <- data[eval(naPat[[2]], data), , drop = FALSE]
+ ## }
+ ## naAct <- eval(mCall$na.action)
+ ## if (!is.null(naAct)) {
+ ## data <- naAct(data)
+ ## }
+ ## subset <- mCall at subset
+ ## if (!is.null(subset)) {
+ ## subset <- eval(asOneSidedFormula(subset)[[2]], data)
+ ## data <- data[subset, ]
+ ## }
+ return(data)
+}
+
+asOneFormula <-
+ ## Constructs a linear formula with all the variables used in a
+ ## list of formulas, except for the names in omit
+ function(..., omit = c(".", "pi"))
+{
+ names <- unique(allVarsRec((list(...))))
+ names <- names[is.na(match(names, omit))]
+ if (length(names)) {
+ eval(parse(text = paste("~", paste(names, collapse = "+")))[[1]])
+ } else NULL
+}
+
+getGroupsFormula <-
+ ## Return the formula(s) for the groups associated with object.
+ ## The result is a one-sided formula unless asList is TRUE in which case
+ ## it is a list of formulas, one for each level.
+ function(object, asList = FALSE, sep = "+")
+ UseMethod("getGroupsFormula")
+
+
+
+getGroupsFormula.default <-
+ ## Return the formula(s) for the groups associated with object.
+ ## The result is a one-sided formula unless asList is TRUE in which case
+ ## it is a list of formulas, one for each level.
+ function(object, asList = FALSE, sep = "/")
+{
+ form <- formula(object)
+ if (!inherits(form, "formula")){
+ stop("\"Form\" argument must be a formula")
+ }
+ form <- form[[length(form)]]
+ if (!((length(form) == 3) && (form[[1]] == as.name("|")))) {
+ ## no conditioning expression
+ return(NULL)
+ }
+ ## val <- list( asOneSidedFormula( form[[ 3 ]] ) )
+ val <- splitFormula(asOneSidedFormula(form[[3]]), sep = sep)
+ names(val) <- unlist(lapply(val, function(el) deparse(el[[2]])))
+# if (!missing(level)) {
+# if (length(level) == 1) {
+# return(val[[level]])
+# } else {
+# val <- val[level]
+# }
+# }
+ if (asList) as.list(val)
+ else as.formula(eval(parse(text = paste("~", paste(names(val),
+ collapse = sep)))))
+}
+
+getGroupsFormula.merMod <- function(object,asList=FALSE, sep="+") {
+ if (asList) {
+ lapply(names(object at flist),asOneSidedFormula)
+ } else {
+ asOneSidedFormula(paste(names(object at flist),collapse=sep))
+ }
+ }
+
+getCovariateFormula <- function (object)
+{
+ form <- formula(object)
+ if (!(inherits(form, "formula"))) {
+ stop("formula(object) must return a formula")
+ }
+ form <- form[[length(form)]]
+ if (length(form) == 3 && form[[1]] == as.name("|")) {
+ form <- form[[2]]
+ }
+ eval(substitute(~form))
+ }
+
+getResponseFormula <-
+ function(object)
+{
+ ## Return the response formula as a one sided formula
+ form <- formula(object)
+ if (!(inherits(form, "formula") && (length(form) == 3))) {
+ stop("\"Form\" must be a two sided formula")
+ }
+ eval(parse(text = paste("~", deparse(form[[2]]))))
+}
+
+##--- needs Trellis/Lattice :
+plot.merMod <-
+ function(x, form = resid(., type = "pearson") ~ fitted(.), abline,
+ id = NULL, idLabels = NULL, idResType = c("pearson", "normalized"),
+ grid, ...)
+ ## Diagnostic plots based on residuals and/or fitted values
+{
+ object <- x
+ if (!inherits(form, "formula")) {
+ stop("\"Form\" must be a formula")
+ }
+ ## constructing data
+ ## can I get away with using object at frame???
+ allV <- all.vars(asOneFormula(form, id, idLabels))
+ allV <- allV[is.na(match(allV,c("T","F","TRUE","FALSE")))]
+ if (length(allV) > 0) {
+ data <- getData(object)
+ if (is.null(data)) { # try to construct data
+ alist <- lapply(as.list(allV), as.name)
+ names(alist) <- allV
+ alist <- c(list(as.name("data.frame")), alist)
+ mode(alist) <- "call"
+ data <- eval(alist, sys.parent(1))
+ } else {
+ if (any(naV <- is.na(match(allV, names(data))))) {
+ stop(paste(allV[naV], "not found in data"))
+ }
+ }
+ } else data <- NULL
+
+ ## this won't do because there may well be variables we want
+ ## that were not in the model call
+
+ ## data <- object at frame
+
+ ## argument list
+ dots <- list(...)
+ if (length(dots) > 0) args <- dots
+ else args <- list()
+ ## appending object to data
+ data <- as.list(c(as.list(data), . = list(object)))
+ ## covariate - must always be present
+ covF <- getCovariateFormula(form)
+ .x <- eval(covF[[2]], data)
+ if (!is.numeric(.x)) {
+ stop("Covariate must be numeric")
+ }
+ argForm <- ~ .x
+ argData <- data.frame(.x = .x, check.names = FALSE)
+ if (is.null(xlab <- attr(.x, "label"))) {
+ xlab <- deparse(covF[[2]])
+ }
+ if (is.null(args$xlab)) args$xlab <- xlab
+
+ ## response - need not be present
+ respF <- getResponseFormula(form)
+ if (!is.null(respF)) {
+ .y <- eval(respF[[2]], data)
+ if (is.null(ylab <- attr(.y, "label"))) {
+ ylab <- deparse(respF[[2]])
+ }
+ argForm <- .y ~ .x
+ argData[, ".y"] <- .y
+ if (is.null(args$ylab)) args$ylab <- ylab
+ }
+
+ ## groups - need not be present
+ grpsF <- getGroupsFormula(form)
+ if (!is.null(grpsF)) {
+ ## ?? FIXME ???
+ gr <- splitFormula(grpsF, sep = "*")
+ for(i in 1:length(gr)) {
+ auxGr <- all.vars(gr[[i]])
+ for(j in auxGr) {
+ argData[[j]] <- eval(as.name(j), data)
+ }
+ }
+ if (length(argForm) == 2)
+ argForm <- eval(parse(text = paste("~ .x |", deparse(grpsF[[2]]))))
+ else argForm <- eval(parse(text = paste(".y ~ .x |", deparse(grpsF[[2]]))))
+ }
+ ## adding to args list
+ args <- c(list(argForm, data = argData), args)
+ if (is.null(args$strip)) {
+ args$strip <- function(...) strip.default(..., style = 1)
+ }
+ if (is.null(args$cex)) args$cex <- par("cex")
+ if (is.null(args$adj)) args$adj <- par("adj")
+
+ if (!is.null(id)) { # identify points in plot
+ idResType <- match.arg(idResType)
+ id <-
+ switch(mode(id),
+ numeric = {
+ if ((id <= 0) || (id >= 1)) {
+ stop("Id must be between 0 and 1")
+ }
+ as.logical(abs(resid(object, type = idResType)) >
+ -qnorm(id / 2))
+ },
+ call = eval(asOneSidedFormula(id)[[2]], data),
+ stop("\"Id\" can only be a formula or numeric.")
+ )
+ if (is.null(idLabels)) {
+ idLabels <- getGroups(object)
+ if (length(idLabels) == 0) idLabels <- 1:object$dims$N
+ idLabels <- as.character(idLabels)
+ } else {
+ if (mode(idLabels) == "call") {
+ idLabels <-
+ as.character(eval(asOneSidedFormula(idLabels)[[2]], data))
+ } else if (is.vector(idLabels)) {
+ if (length(idLabels <- unlist(idLabels)) != length(id)) {
+ stop("\"IdLabels\" of incorrect length")
+ }
+ idLabels <- as.character(idLabels)
+ } else {
+ stop("\"IdLabels\" can only be a formula or a vector")
+ }
+ }
+ }
+
+ ## defining abline, if needed
+ if (missing(abline)) {
+ if (missing(form)) { # r ~ f
+ abline <- c(0, 0)
+ } else {
+ abline <- NULL
+ }
+ }
+
+ #assign("id", id , where = 1)
+ #assign("idLabels", idLabels, where = 1)
+ #assign("abl", abline, where = 1)
+ assign("abl", abline)
+
+ ## defining the type of plot
+ if (length(argForm) == 3) {
+ if (is.numeric(.y)) { # xyplot
+ plotFun <- "xyplot"
+ if (is.null(args$panel)) {
+ args <- c(args,
+ panel = list(function(x, y, subscripts, ...)
+ {
+ x <- as.numeric(x)
+ y <- as.numeric(y)
+ dots <- list(...)
+ if (grid) panel.grid()
+ panel.xyplot(x, y, ...)
+ if (any(ids <- id[subscripts])){
+ ltext(x[ids], y[ids], idLabels[subscripts][ids],
+ cex = dots$cex, adj = dots$adj)
+ }
+ if (!is.null(abl)) {
+ if (length(abl) == 2) panel.abline(a = abl, ...) else panel.abline(h = abl, ...)
+ }
+ }))
+ }
+ } else { # assume factor or character
+ plotFun <- "bwplot"
+ if (is.null(args$panel)) {
+ args <- c(args,
+ panel = list(function(x, y, ...)
+ {
+ if (grid) panel.grid()
+ panel.bwplot(x, y, ...)
+ if (!is.null(abl)) {
+ panel.abline(v = abl[1], ...)
+ }
+ }))
+ }
+ }
+ } else {
+ plotFun <- "histogram"
+ if (is.null(args$panel)) {
+ args <- c(args,
+ panel = list(function(x, ...)
+ {
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/lme4 -r 1644
More information about the Lme4-commits
mailing list