[Lme4-commits] r1688 - pkg/lme4/R

[noreply at r-forge.r-project.org]

Author: bbolker
Date: 2012-03-22 19:41:39 +0100 (Thu, 22 Mar 2012)
New Revision: 1688

Modified:
   pkg/lme4/R/lmer.R
   pkg/lme4/R/profile.R
   pkg/lme4/R/utilities.R
Log:

   made refitML use optimizer argument/optwrap
   tweaked optwrap (less dependence on rho)
   added comments about first steps in GLM fitting procedure
   first steps toward more robust profile (induces warnings
     about code/doc mismatch in profile ...)



Modified: pkg/lme4/R/lmer.R
===================================================================
--- pkg/lme4/R/lmer.R	2012-03-22 18:30:20 UTC (rev 1687)
+++ pkg/lme4/R/lmer.R	2012-03-22 18:41:39 UTC (rev 1688)
@@ -1115,7 +1115,9 @@
 }
 
 ##' @S3method refitML merMod
-refitML.merMod <- function (x, ...) {
+refitML.merMod <- function (x, optimizer="bobyqa", ...) {
+    ## FIXME: optimizer is set to 'bobyqa' for back-compatibility, but that's not
+    ##  consistent with lmer (default NM).  Should be based on internally stored 'optimizer' value
     if (!isREML(x)) return(x)
     stopifnot(is(rr <- x at resp, "lmerResp"))
     rho <- new.env(parent=parent.env(environment()))
@@ -1124,7 +1126,8 @@
     rho$pp <- new(class(xpp), X=xpp$X, Zt=xpp$Zt, Lambdat=xpp$Lambdat,
                   Lind=xpp$Lind, theta=xpp$theta, n=nrow(xpp$X))
     devfun <- mkdevfun(rho, 0L)
-    opt <- bobyqa(x at theta, devfun, x at lower)
+    opt <- optwrap(optimizer, devfun, x at theta, lower=x at lower)
+    ##  opt <- bobyqa(x at theta, devfun, x at lower)
     n <- length(rr$y)
     pp <- rho$pp
     p <- ncol(pp$X)
@@ -1948,11 +1951,12 @@
 }
 
 optwrap <- function(optimizer, fn, par, lower=-Inf, upper=Inf,
-                    control=list(), rho, adj=FALSE, verbose=0L) {
-  ## control and rho must be specified if adj==TRUE;
-  ##  otherwise this is a fairly simple wrapper
-  optfun <- getOptfun(optimizer)
-  ## special-purpose control parameter tweaks: only for second round in nlmer, glmer
+                    control=list(), rho=NULL, adj=FALSE, verbose=0L) {
+    ## control and rho must be specified if adj==TRUE;
+    ##  otherwise this is a fairly simple wrapper
+    ## FIXME: would like to remove rho-dependence (used to pass back modified control settings) ?
+    optfun <- getOptfun(optimizer)
+    ## special-purpose control parameter tweaks: only for second round in nlmer, glmer
 
   lower <- rep(lower, length.out=length(par))
   upper <- rep(upper, length.out=length(par))
@@ -1962,18 +1966,18 @@
            bobyqa = {
              if(!is.numeric(control$rhobeg)) control$rhobeg <- 0.0002
              if(!is.numeric(control$rhoend)) control$rhoend <- 2e-7
-             rho$control <- control
+             if (!is.null(rho)) 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())))
+                 xst <- c(rep.int(0.1, length(environment(fn)$pp$theta)),  ## theta parameters
+                          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, upper=upper, control=control)
+             if (!is.null(rho)) rho$control <- control
+         })
+    arglist <- list(fn=fn, par=par, lower=lower, upper=upper, 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?

Modified: pkg/lme4/R/profile.R
===================================================================
--- pkg/lme4/R/profile.R	2012-03-22 18:30:20 UTC (rev 1687)
+++ pkg/lme4/R/profile.R	2012-03-22 18:41:39 UTC (rev 1688)
@@ -46,6 +46,7 @@
 profile.merMod <- function(fitted, which=1:nptot, alphamax = 0.01, maxpts = 100, delta = cutoff/8,
                            ##  tr = 0,  ## FIXME:  remove if not doing anything ...
                            verbose=0, devtol=1e-9,
+                           maxmult = 10,
                            startmethod = "prev",
                            optimizer="bobyqa", ...) {
 
@@ -96,7 +97,7 @@
     ## (absstep) and the values of zeta and column cc for rows
     ## r-1 and r.  The parameter may not be below lower (or above upper)
     nextpar <- function(mat, cc, r, absstep,
-                        lower = -Inf, upper = Inf, minstep=1e-6, maxmult=10.0) {
+                        lower = -Inf, upper = Inf, minstep=1e-6) {
         rows <- r - (1:0)         # previous two row numbers
         pvals <- mat[rows, cc]
         zeta <- mat[rows, ".zeta"]

Modified: pkg/lme4/R/utilities.R
===================================================================
--- pkg/lme4/R/utilities.R	2012-03-22 18:30:20 UTC (rev 1687)
+++ pkg/lme4/R/utilities.R	2012-03-22 18:41:39 UTC (rev 1688)
@@ -134,6 +134,7 @@
 ##' @family utilities
 ##' @export
 mkRespMod <- function(fr, REML=NULL, family = NULL, nlenv = NULL, nlmod = NULL) {
+    ## FIXME: may need to add X, or pass it somehow, if we want to use glm.fit
     y <- model.response(fr)
     if(length(dim(y)) == 1) {
 	## avoid problems with 1D arrays, but keep names
@@ -181,6 +182,9 @@
     }
     stopifnot(inherits(family, "family"))
                               # need weights for initialize evaluation
+    ## FIXME: family$initialize may not be good enough, need
+    ##        glm.fit?
+    ## do.call(glm.fit,c(as.list(rho),list(x=X,family=family)))
     rho$nobs <- n
     eval(family$initialize, rho)
     family$initialize <- NULL     # remove clutter from str output