[Lme4-commits] r1607 - in pkg/lme4Eigen: R tests

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

Author: bbolker
Date: 2012-02-21 17:12:33 +0100 (Tue, 21 Feb 2012)
New Revision: 1607

Modified:
   pkg/lme4Eigen/R/lmer.R
   pkg/lme4Eigen/tests/refit.R
Log:

  further refit/simulate cleanup
  add isGLMM() method
  add Poisson tests to refit.R



Modified: pkg/lme4Eigen/R/lmer.R
===================================================================
--- pkg/lme4Eigen/R/lmer.R	2012-02-21 16:04:28 UTC (rev 1606)
+++ pkg/lme4Eigen/R/lmer.R	2012-02-21 16:12:33 UTC (rev 1607)
@@ -1165,14 +1165,19 @@
 ##' @S3method refit merMod
 refit.merMod <- function(object, newresp, ...) {
     rr        <- object at resp$copy()
-    ## FIXME: want this to work for binomial (two-column) data?
-    ##  or tell people they have to work with prop/weights formulation?
-    if (is.matrix(newresp) && ncol(newresp)==2 && rr$family$family=="binomial") {
-      ntot <- rowSums(newresp)
-      ## FIXME: maybe unnecessary?
-      if (any(ntot==0)) stop("totals must be >0 ")
-      newresp <- newresp[,1]/ntot
-      rr$setWeights(ntot)
+
+    if (isGLMM(object) && rr$family$family=="binomial") {
+      if (is.matrix(newresp) && ncol(newresp)==2) {
+        ntot <- rowSums(newresp)
+        ## FIXME: test what happens for (0,0) columns
+        newresp <- newresp[,1]/ntot
+        rr$setWeights(ntot)
+      }
+      if (is.factor(newresp)) {
+        ## FIXME: would be better to do this consistently with
+        ## whatever machinery is used in glm/glm.fit/glmer ... ??
+        newresp <- as.numeric(newresp)-1
+      }
     }
     stopifnot(length(newresp <- as.numeric(as.vector(newresp))) == length(rr$y))
     rr$setResp(newresp)
@@ -1180,10 +1185,13 @@
     dc        <- object at devcomp
     nAGQ      <- dc$dims["nAGQ"]
     nth       <- dc$dims["nth"]
-    ff <- mkdevfun(list2env(list(pp=pp, resp=rr, u0=pp$u0, beta0=pp$beta0, verbose=0L,
-                                 tolPwrss=dc$cmp["tolPwrss"], dpars=seq_len(nth),
-                                 nAGQ=unname(nAGQ))),
-                        nAGQ=nAGQ)
+    devlist <- list(pp=pp, resp=rr, u0=pp$u0, beta0=pp$beta0, verbose=0L,
+                    dpars=seq_len(nth))
+    if (isGLMM(object)) {
+      devlist <- c(devlist,tolPwrss=unname(dc$cmp["tolPwrss"]),
+                   nAGQ=unname(nAGQ))
+    }
+    ff <- mkdevfun(list2env(devlist),nAGQ=nAGQ)
     xst       <- rep.int(0.1, nth)
     x0        <- pp$theta
     lower     <- object at lower
@@ -1197,7 +1205,7 @@
 ### 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)
     mkMerMod(environment(ff), opt, list(flist=object at flist, cnms=object at cnms, Gp=object at Gp,
-                                        lower=lower), object at frame, getCall(object))
+                                        lower=object at lower), object at frame, getCall(object))
 }    
 
 ##' @S3method refitML merMod
@@ -1981,3 +1989,8 @@
 weights.merMod <- function(object, ...) {
   object at resp$weights
 }
+
+isGLMM <- function(object) {
+  as.logical(object at devcomp$dims["GLMM"])
+}
+             

Modified: pkg/lme4Eigen/tests/refit.R
===================================================================
--- pkg/lme4Eigen/tests/refit.R	2012-02-21 16:04:28 UTC (rev 1606)
+++ pkg/lme4Eigen/tests/refit.R	2012-02-21 16:12:33 UTC (rev 1607)
@@ -60,6 +60,18 @@
 getinfo(gm3R)
 getinfo(gm3S)
 
+data(Mmmec,package="mlmRev")
+gm4 <- lmer(deaths ~ uvb + (1|region), data=Mmmec,
+            family=poisson,
+            offset = log(expected))
+gm4R <- refit(gm4,Mmmec$deaths)
+gm4S <- refit(gm4,simulate(gm4)[[1]])
+getinfo(gm4)
+getinfo(gm4R)
+getinfo(gm4S)
+
+stopifnot(all.equal(gm4,gm4R,tol=5e-5))
+
 checkcomp <- function(x,y) {
   for (i in slotNames(x)) {
     cat(i,"\n")