[R-sig-ME] simulate glmm with family=binomial
Austin Frank
austin.frank at gmail.com
Tue Jun 17 06:33:31 CEST 2008
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
More information about the R-sig-mixed-models
mailing list