[R-sig-ME] simulate() fails for poisson lmer fits
William Valdar
valdar at well.ox.ac.uk
Thu Aug 2 18:50:53 CEST 2007
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:
>
> 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:
>> 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
>>
>>
>
> =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
> 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
>
>
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
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
More information about the R-sig-mixed-models
mailing list