Skip to content

simulate() fails for poisson lmer fits

5 messages · William Valdar, Martin Maechler

Original
William Valdar
# 
Hello,

I wish to simulate() from a fitted Poisson GLMM. Both lmer() and lmer2() 
from lme4 (version info at the bottom) fail, apparently due to bugs. 
Here's a test 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)
# Error: inherits(object, "lmer") is not TRUE

fit <- lmer2(count ~ weight + (1|batch),  family=poisson, data=counts)
simulate(fit)
# CHOLMOD error: X and/or Y have wrong dimensions
# Error in crossprod(object at ZXyt, c(unlist(lapply(seq_along(re), 
# function(k) (t(cholL[[k]]) %*%  :
#        Cholmod error 'X and/or Y have wrong dimensions' at 
# file:../MatrixOps/cholmod_sdmult.c, line 90

Is there a quick fix for either of these two? Otherwise, is there an 
alternative (I've checked objects produced by nlme, glmmPQL, GLMMGibbs 
with no luck)? I am using lme4 0.99875-6 version R 2.5.1 on Windows XP.

Many thanks,

Will

=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
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
William Valdar
# 
I missed out a line in the code below. The line that throws the error is 
not the fitting but the subsequent simulation, ie:

fit <- lmer(count ~ weight + (1|batch),  family=poisson, data=counts)
simulate(fit)
# Error: inherits(object, "lmer") is not TRUE
On Wed, 1 Aug 2007, William Valdar wrote:
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
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
Martin Maechler
# 
Hi,
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
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:
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
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
William Valdar
# 
I spoke too soon. The simulate() runs but produces a real-valued response 
rather than an integral response. Here's a modification of the code you 
sent for gaussian and poisson responses:

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
     lpred <- lpred + drop(object at X %*% fixef(object))
     n <- prod(dim(lpred))

     response <- NULL
     if (is.null(attr(object, "family")))
     {
          ## and per-observation noise term
          response <- as.data.frame(lpred + rnorm(n, sd = sc))
     }
     else
     {
          fam <- attr(object, "family")
          if ("poisson"==fam$family)
          {
              response <- as.data.frame(matrix(byrow = FALSE, ncol=nsim,
                      rpois(n, fam$linkinv(lpred))))
          }
          else
          {
              stop("simulate() not yet implemented for", fam$family,
                      "glmms\n")
          }
     }

     ## save the seed
     attr(response, "seed") <- RNGstate
     response
     })
On Thu, 2 Aug 2007, William Valdar wrote:
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
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