[R-sig-ME] simulate() fails for poisson lmer fits

[William Valdar]

That works perfectly, thanks very much!

It's also considerably faster than the hack fix I wrote in the early hours 
of this morning (included below only for interest and as a warning to 
others).

Will

---hack---

rnorm.factor <- function(g, sd=1, mean=0)
# generate rnorm for a simple "batch" effect
{
     g   <- as.factor(g)[drop=TRUE]
     x.g <- rnorm( nlevels(g), sd=sd, mean=mean )
     return (x.g[g])
}

simulate.lmer.hack <- function(fit, nsim=1, seed=NULL, ...)
# simulate from fitted lmer obj with simple random intercepts
{
     if (!is.null(seed)) set.seed(seed)

     n <- length(fitted(fit))
     fixed.linpreds <- model.matrix(attr(fit, "terms"),
             data=attr(fit, "frame")) %*% fixef(fit)

     ymat <- matrix(NA, nrow=n, ncol=nsim)
     for (s in 1:nsim)
     {
         random.linpreds <- rep(0, length(n))
         vc <- VarCorr(fit)
         for (i in 1:length(vc))
         {
             random.sd   <- sqrt(vc[[i]][1])
             random.data <- attr(fit, "flist")[[i]]
             random.linpreds <- random.linpreds +
                     rnorm.factor(random.data, sd=random.sd, mean=0)
         }
         # make canonical param
         theta <- fixed.linpreds + random.linpreds
         if (is.null(attr(fit, "family")))
         {
             scale <- attr(vc, "sc")
             ymat[,s] <- rnorm(n, theta, scale)
         }
         else if ("poisson" == attr(fit, "family")$family)
         {
             ymat[,s] <- rpois(n, attr(fit, "family")$linkinv(theta))
         }
     }
     as.data.frame(ymat)
}

counts <- data.frame(
         count  = c(8, 3, 0, 3, 0, 9, 0, 11, 4, 7, 0, 0, 0, 4, 3, 6, 3,
                 15, 11, 9),
         batch = c(3.1, 3.1, 3.1, 3.3, 3.3, 3.3, 3.2, 3.12, 3.8, 3.11,
                 3.4, 3.4, 3.4, 3.4, 3.5, 3.5, 3.5, 3.5, 3.6, 3.6),
         weight = c(324.4, 372.5, 352.7, 379.6, 388.1, 431, 448.4, 377.3,
                 376.5, 358.4, 356, 351.4, 350.8, 332.1, 334.5, 392, 370.5,
                 409.7, 375, 318.5))

fit.poisson <- lmer(count ~ weight + (1|batch),  family=poisson, 
data=counts)
simulate.lmer.hack(fit.poisson)


On Thu, 2 Aug 2007, Martin Maechler wrote:
> Hi,
>>>>>> "WV" == William Valdar <valdar at well.ox.ac.uk>
>>>>>>     on Wed, 1 Aug 2007 18:49:58 +0100 (BST) writes:
>
>    WV> Hello, I wish to simulate() from a fitted Poisson
>    WV> GLMM. Both lmer() and lmer2() from lme4 (version info at
>    WV> the bottom) fail, apparently due to bugs.  Here's a test
>    WV> case:
>
>    counts <- data.frame(
> 	     count  = c(8, 3, 0, 3, 0, 9, 0, 11, 4, 7, 0, 0, 0, 4, 3, 6, 3,
> 		     15, 11, 9),
> 	     batch = c(3.1, 3.1, 3.1, 3.3, 3.3, 3.3, 3.2, 3.12, 3.8, 3.11,
> 		     3.4, 3.4, 3.4, 3.4, 3.5, 3.5, 3.5, 3.5, 3.6, 3.6),
> 	     weight = c(324.4, 372.5, 352.7, 379.6, 388.1, 431, 448.4, 377.3,
> 		     376.5, 358.4, 356, 351.4, 350.8, 332.1, 334.5, 392, 370.5,
> 		     409.7, 375, 318.5))
>
>    fit <- lmer(count ~ weight + (1|batch),  family=poisson, data=counts)
>
> MM, that works fine, but
>
>    simulate(fit)
>
> gives what you say:
>
>  WV>    # Error: inherits(object, "lmer") is not TRUE
>
> and that's indeed a bug in the simulate method:
>
> Using  inherits(object, "lmer")  is wrong and should be
> replaced by  is(object, "lmer") .
>
> Instead of waiting for the next version of lme4,
> I think the following should work for you:
>
>
> setMethod("simulate", signature(object = "mer"),
> 	  function(object, nsim = 1, seed = NULL, ...)
>      {
> 	  if(!exists(".Random.seed", envir = .GlobalEnv))
> 	      runif(1)		     # initialize the RNG if necessary
> 	  if(is.null(seed))
> 	      RNGstate <- .Random.seed
> 	  else {
> 	      R.seed <- .Random.seed
> 	      set.seed(seed)
> 	      RNGstate <- structure(seed, kind = as.list(RNGkind()))
> 	      on.exit(assign(".Random.seed", R.seed, envir = .GlobalEnv))
> 	  }
>
> 	  stopifnot((nsim <- as.integer(nsim[1])) > 0, is(object, "lmer"))
> 	  ## similate the linear predictors
> 	  lpred <- .Call(lme4:::mer_simulate, object, nsim)
> 	  sc <- abs(object at devComp[8])
>
> 	  ## add fixed-effects contribution and per-observation noise term
> 	  lpred <- as.data.frame(lpred + drop(object at X %*% fixef(object)) +
> 				 rnorm(prod(dim(lpred)), sd = sc))
> 	  ## save the seed
> 	  attr(lpred, "seed") <- RNGstate
> 	  lpred
>      })
>
>
> at least it fixes the problem for me.
>
>
>    WV> fit <- lmer2(count ~ weight + (1|batch),  family=poisson, data=counts)
>    WV> simulate(fit)
>    WV> # CHOLMOD error: X and/or Y have wrong dimensions
>    WV> # Error in crossprod(object at ZXyt, c(unlist(lapply(seq_along(re),
>    WV> # function(k) (t(cholL[[k]]) %*%  :
>    WV> #        Cholmod error 'X and/or Y have wrong dimensions' at
>    WV> # file:../MatrixOps/cholmod_sdmult.c, line 90
>
> I can confirm that error and I agree that it is a bug,
> well at least in the error message :-)
>
> BTW, also     summary(fit)   gives an error.
> {which is caused by  vcov(fit)  getting 0 x 0 matrix }
>
> But then,  lmer2() is in beta testing,
> so thanks a lot for your nicely reproducible example.
>
>
> Note that your data set has only one observation for some batch
> levels which may be deemed to give problems,
> but in fact that's not the problem here at all.
>
> Regards,
> Martin Maechler, ETH Zurich
>
>
>    WV> Is there a quick fix for either of these two? Otherwise, is there an
>    WV> alternative (I've checked objects produced by nlme, glmmPQL, GLMMGibbs
>    WV> with no luck)? I am using lme4 0.99875-6 version R 2.5.1 on Windows XP.
>
>    WV> Many thanks,
>
>    WV> Will
>
>    WV> =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
>    WV> Dr William Valdar               ++44 (0)1865 287 589
>    WV> Wellcome Trust Centre           valdar at well.ox.ac.uk
>    WV> for Human Genetics, Oxford      www.well.ox.ac.uk/~valdar
>
>

=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
Dr William Valdar               ++44 (0)1865 287 589
Wellcome Trust Centre           valdar at well.ox.ac.uk
for Human Genetics, Oxford      www.well.ox.ac.uk/~valdar