[R-sig-ME] simulate glmm with family=binomial

Austin Frank austin.frank at gmail.com
Tue Jun 17 06:33:31 CEST 2008


In lme4 version 0.99875-9, simulate() is implemented for linear mixed
models and generalized linear mixed models from the Poisson family.  I'm
hoping to simulate a glmm from the binomial family.

I looked at the code used to implement the Poisson simulations, and
naively tried to port it to handle the binomial (and only the binomial)
case.  The result is:

--8<---------------cut here---------------start------------->8---
my.simulate <- function (object, nsim = 1, seed = NULL, ...) 
  if (!exists(".Random.seed", envir = .GlobalEnv)) 
  if (is.null(seed)) 
    RNGstate <- .Random.seed
  else {
    R.seed <- .Random.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"))
  lpred <- .Call(lme4:::mer_simulate, object, nsim)
  sc <- abs(object at devComp[8])
  lpred <- lpred + drop(object at X %*% fixef(object))
  n <- prod(dim(lpred))
  fam <- attr(object, "family")
  if (fam$family == "binomial") {
    response <- as.data.frame(matrix(byrow = FALSE, ncol = nsim, 
                                     rbinom(n, size=1, plogis(lpred))))
    attr(response, "seed") <- RNGstate
--8<---------------cut here---------------end--------------->8---

The part I changed is:

  response <- as.data.frame(matrix(byrow = FALSE, ncol = nsim, 
                                   rbinom(n, size=1, plogis(lpred))))

which is based on (ripped off from) the code for Poisson simulations:

  response <- as.data.frame(matrix(byrow = FALSE, ncol = nsim, 
                                   rpois(n, fam$linkinv(lpred))))

This works, in the sense that I get the correct number of samples, each
with a value of 0 or 1.  But I have to assume that there's a problem
here-- if this were all it took, it would be included already!
So, three questions:
1.  Is there something wrong with this solution?
2.  Assuming so, is there an obvious fix to the problem?
3.  In the second argument to rbinom, I've specified that we're basing
    this on a single trial.  Am I right in thinking that size=1 is the
    right argument here?

Thanks in advance for any help!
Austin Frank
GPG Public Key (D7398C2F): http://aufrank.net/personal.asc

More information about the R-sig-mixed-models mailing list