[Lme4-commits] r1774 - in pkg/lme4: R tests

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

Author: bbolker
Date: 2012-06-24 10:19:14 +0200 (Sun, 24 Jun 2012)
New Revision: 1774

Modified:
   pkg/lme4/R/bootMer.R
   pkg/lme4/R/lmer.R
   pkg/lme4/tests/simulate.R
Log:

  slightly extend error message in bootMer
  fix use.u in simulate (remaining uncertainty/ambiguity in use.u definition!)



Modified: pkg/lme4/R/bootMer.R
===================================================================
--- pkg/lme4/R/bootMer.R	2012-06-23 11:27:12 UTC (rev 1773)
+++ pkg/lme4/R/bootMer.R	2012-06-24 08:19:14 UTC (rev 1774)
@@ -120,7 +120,7 @@
             ss <- replicate(nsim,fitted(x)+sample(residuals(x,"response")),
                             simplify=FALSE)
         } else {
-            stop("not implemented")
+            stop("semiparametric bootstrapping with use.u=FALSE not yet implemented")
         }
     }
     t.star <- matrix(t0, nrow = length(t0), ncol = nsim)

Modified: pkg/lme4/R/lmer.R
===================================================================
--- pkg/lme4/R/lmer.R	2012-06-23 11:27:12 UTC (rev 1773)
+++ pkg/lme4/R/lmer.R	2012-06-24 08:19:14 UTC (rev 1774)
@@ -1310,8 +1310,9 @@
 ##' @param nsim positive integer scalar - the number of responses to simulate
 ##' @param seed an optional seed to be used in \code{set.seed} immediately
 ##'     before the simulation so as to generate a reproducible sample.
-##' @param use.u (logical) generate new random-effects values (FALSE) or
-##'     generate a simulation condition on the current random-effects estimates (TRUE)?
+##' @param use.u (logical) generate a simulation conditional on the current
+##' random-effects estimates (TRUE) or generate new random-effects values (FALSE) [FIXME!?]
+
 ##' @param ... optional additional arguments, none are used at present
 ##' @examples
 ##' ## test whether fitted models are consistent with the
@@ -1342,14 +1343,14 @@
     if (length(offset <- getME(object,"offset"))>0) {
       etasim.fix <- etasim.fix+offset
     }
+    U <- getME(object, "Z") %*% getME(object, "Lambda")
+    u <- if (use.u) {
+        rep(getME(object, "u"), nsim)
+    } else {
+        rnorm(ncol(U)*nsim)
+    }
     etasim.reff <- ## UNSCALED random-effects contribution:
-      if(use.u) {
-        getME(object, "u")
-      } else {
-        U <- getME(object, "Z") %*% getME(object, "Lambda")
-        q <- ncol(U)
-        as(U %*% matrix(rnorm(q * nsim), ncol = nsim), "matrix")
-      }
+        as(U %*% matrix(u, ncol = nsim), "matrix")
     if (is(object at resp,"lmerResp")) {
       ## result will be matrix  n x nsim :
       val <- etasim.fix + sigma * (etasim.reff +

Modified: pkg/lme4/tests/simulate.R
===================================================================
--- pkg/lme4/tests/simulate.R	2012-06-23 11:27:12 UTC (rev 1773)
+++ pkg/lme4/tests/simulate.R	2012-06-24 08:19:14 UTC (rev 1774)
@@ -1,5 +1,9 @@
 library(lme4)
 
+fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
+s1 <- simulate(fm1,seed=101)[[1]]
+s2 <- simulate(fm1,seed=101,use.u=TRUE)
+
 ## binomial (2-column and prob/weights)
 gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
               data = cbpp, family = binomial)
@@ -9,11 +13,21 @@
 s1 <- simulate(gm1,seed=101)[[1]]
 s2 <- simulate(gm2,seed=101)[[1]]
 stopifnot(all.equal(s1[,1]/rowSums(s1),s2))
+s3 <- simulate(gm1,seed=101,use.u=TRUE)
 
+## test explicitly stated link function
 gm3 <- glmer(cbind(incidence, size - incidence) ~ period +
              (1 | herd), data = cbpp, family = binomial(link="logit"))
-s3 <- simulate(gm3,seed=101)[[1]]
-stopifnot(all.equal(s3,s1))
+s4 <- simulate(gm3,seed=101)[[1]]
+stopifnot(all.equal(s1,s4))
+
+cbpp$obs <- factor(seq(nrow(cbpp)))
+gm4 <- glmer(cbind(incidence, size - incidence) ~ period +
+             (1 | herd) + (1|obs), data = cbpp, family = binomial)
+
+s5 <- simulate(gm4,seed=101)[[1]]
+s6 <- simulate(gm4,seed=101,use.u=TRUE)[[1]]
+
 ## Bernoulli
 ## works, but too slow
 if (FALSE) {