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

[Martin Maechler]

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