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

[Austin Frank]

Hello!

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)) 
    runif(1)
  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"))
  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") {
    summary(lpred)
    response <- as.data.frame(matrix(byrow = FALSE, ncol = nsim, 
                                     rbinom(n, size=1, plogis(lpred))))
    attr(response, "seed") <- RNGstate
    return(response)
  }
}
--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!
/au
    
-- 
Austin Frank
http://aufrank.net
GPG Public Key (D7398C2F): http://aufrank.net/personal.asc