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

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

Author: bbolker
Date: 2012-07-08 21:12:21 +0200 (Sun, 08 Jul 2012)
New Revision: 1785

Modified:
   pkg/lme4/R/lmer.R
   pkg/lme4/R/predict.R
Log:

  added boundary-value tests to lmer (*not* glmer etc. ... modularize!)
  added na.action *stub* to predict.R



Modified: pkg/lme4/R/lmer.R
===================================================================
--- pkg/lme4/R/lmer.R	2012-06-28 10:10:02 UTC (rev 1784)
+++ pkg/lme4/R/lmer.R	2012-07-08 19:12:21 UTC (rev 1785)
@@ -82,7 +82,7 @@
                  optimizer="Nelder_Mead", ...)
 {
     verbose <- as.integer(verbose)
-    restart <- FALSE
+    restart <- TRUE ## FIXME; set default elsewhere?
     if (!is.null(control$restart)) {
         restart <- control$restart
         control$restart <- NULL
@@ -168,11 +168,35 @@
                    devfun, rho$pp$theta, lower=reTrms$lower, control=control,
                    adj=FALSE, verbose=verbose)
 
-    if (restart && any(rho$pp$theta==0)) {
-        if (verbose) message("some theta parameters on the boundary, restarting")
-        opt <- optwrap(optimizer,
-                       devfun, rho$pp$theta, lower=reTrms$lower, control=control,
-                       adj=FALSE, verbose=verbose)
+    if (restart) {
+        ## FIXME: should we be looking at rho$pp$theta or opt$par
+        ##  at this point???  in koller example (for getData(13)) we have
+        ##   rho$pp$theta=0, opt$par=0.08
+        if (length(bvals <- which(rho$pp$theta==reTrms$lower))>0) {
+            ## *don't* use numDeriv -- cruder but fewer dependencies, no worries
+            ##  about keeping to the interior of the allowed space
+            theta0 <- new("numeric",rho$pp$theta) ## 'deep' copy ...
+            d0 <- devfun(theta0)
+            btol <- 1e-5  ## FIXME: make user-settable?
+            ## FIXME: opt$fval is wrong
+            bgrad <- sapply(bvals,
+                            function(i) {
+                                bndval <- reTrms$lower[i]
+                                theta <- theta0
+                                theta[i] <- bndval+btol
+                                (devfun(theta)-d0)/btol
+                            })
+            ## what do I need to do to reset rho$pp$theta to original value???
+            devfun(theta0) ## reset rho$pp$theta after tests
+            ## FIXME: allow user to specify ALWAYS restart if on boundary?
+            if (any(bgrad<0)) {
+                if (verbose) message("some theta parameters on the boundary, restarting")
+                opt <- optwrap(optimizer,
+                               devfun, opt$par,
+                               lower=reTrms$lower, control=control,
+                               adj=FALSE, verbose=verbose)
+            }
+        }
     }
         
     mkMerMod(environment(devfun), opt, reTrms, fr, mc)

Modified: pkg/lme4/R/predict.R
===================================================================
--- pkg/lme4/R/predict.R	2012-06-28 10:10:02 UTC (rev 1784)
+++ pkg/lme4/R/predict.R	2012-07-08 19:12:21 UTC (rev 1785)
@@ -12,6 +12,7 @@
 ##'    detected in \code{newdata} will trigger an error; if TRUE, then
 ##'    the prediction will use the unconditional (population-level)
 ##'    values for data with previously unobserved levels
+##' @param na.action function determining what should be done with missing values in \code{newdata}.  Currently a stub: the default will be to predict ‘NA’. See \code{\link{na.pass}}.
 ##' @param ... optional additional parameters.  None are used at present.
 ##' @return a numeric vector of predicted values
 ##' @note explain why there is no option for computing standard errors of predictions!
@@ -27,8 +28,14 @@
 ##' @export
 predict.merMod <- function(object, newdata=NULL, REform=NULL,
                            terms=NULL, type=c("link","response"),
-                           allow.new.levels=FALSE, ...) {
+                           allow.new.levels=FALSE, na.action=na.pass, ...) {
     ## FIXME: appropriate names for result vector?
+
+    if (length(list(...)>0)) warning("unused arguments ignored")
+    if (!missing(na.action)) {
+        ## FIXME
+        stop("na.action is not yet implemented.")
+    }
     type <- match.arg(type)
     if (!is.null(terms)) stop("terms functionality for predict not yet implemented")
     X_orig <- getME(object, "X")
@@ -45,7 +52,8 @@
             RHS <- formula(object,fixed.only=TRUE)[-2]
             Terms <- terms(object,fixed.only=TRUE)
             X <- model.matrix(RHS, model.frame(delete.response(Terms), newdata),
-                              contrasts.arg=attr(X_orig,"contrasts"))
+                              contrasts.arg=attr(X_orig,"contrasts"),
+                              na.action=attr(object at frame,"na.action"))
         }
         pred <- drop(X %*% fixef(object))
         ## modified from predict.glm ...