proposed simulate.glm method

I have found the "simulate" method (incorporated
in some packages) very handy. As far as I can tell the
only class for which simulate is actually implemented
in base R is lm ... this is actually a little dangerous
for a naive user who might be tempted to try
simulate(X) where X is a glm fit instead, because
it defaults to simulate.lm (since glm inherits from
the lm class), and the answers make no sense ...

Here is my simulate.glm(), which is modeled on
simulate.lm .  It implements simulation for poisson
and binomial (binary or non-binary) models, should
be easy to implement others if that seems necessary.

  I hereby request comments and suggest that it wouldn't
hurt to incorporate it into base R ...  (I will write
docs for it if necessary, perhaps by modifying ?simulate --
there is no specific documentation for simulate.lm)

  cheers
    Ben Bolker

simulate.glm <- function (object, nsim = 1, seed = NULL, ...)
{
  ## RNG stuff copied from simulate.lm
  if (!exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE))
    runif(1)
  if (is.null(seed))
    RNGstate <- get(".Random.seed", envir = .GlobalEnv)
  else {
    R.seed <- get(".Random.seed", envir = .GlobalEnv)
    set.seed(seed)
    RNGstate <- structure(seed, kind = as.list(RNGkind()))
    on.exit(assign(".Random.seed", R.seed, envir = .GlobalEnv))
  }
  ## get probabilities/intensities
  pred <- matrix(rep(predict(object,type="response"),nsim),ncol=nsim)
  ntot <- length(pred)
  if (object$family$family=="binomial") {
    resp <- object$model[[1]]
    size <- if (is.matrix(resp)) rowSums(resp) else 1
  }
  val <- switch(object$family$family,
                poisson=rpois(ntot,pred),
                binomial=rbinom(ntot,prob=pred,size=size),
                stop("family ",object$family$family," not implemented"))
  ans <- as.data.frame(matrix(val,ncol=nsim))
  attr(ans, "seed") <- RNGstate
  ans
}

if (FALSE) {
  ## examples: modified from ?simulate
  x <- 1:10
  n <- 10
  y <- rbinom(length(x),prob=plogis((x-5)/2),size=n)
  y2 <- c("a","b")[1+rbinom(length(x),prob=plogis((x-5)/2),size=1)]
  mod1 <- glm(cbind(y,n-y) ~ x,family=binomial)
  mod2 <- glm(factor(y2) ~ x,family=binomial)
  S1 <- simulate(mod1, nsim = 4)
  S1B <- simulate(mod2, nsim = 4)
  ## repeat the simulation:
  .Random.seed <- attr(S1, "seed")
  identical(S1, simulate(mod1, nsim = 4))

  S2 <- simulate(mod1, nsim = 200, seed = 101)
  rowMeans(S2)/10 # after correcting for binomial sample size, should be
about
  fitted(mod1)

  plot(rowMeans(S2)/10)
  lines(fitted(mod1))

  ## repeat identically:
  (sseed <- attr(S2, "seed")) # seed; RNGkind as attribute
  stopifnot(identical(S2, simulate(mod1, nsim = 200, seed = sseed)))
}