[Rd] proposed simulate.glm method

From: Ben Bolker <bolker_at_ufl.edu>
Date: Thu, 12 Feb 2009 11:29:14 -0500

  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))) }

-- 
Ben Bolker
Associate professor, Biology Dep't, Univ. of Florida
bolker_at_ufl.edu / www.zoology.ufl.edu/bolker
GPG key: www.zoology.ufl.edu/bolker/benbolker-publickey.asc


______________________________________________ R-devel_at_r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel

Received on Thu 12 Feb 2009 - 15:50:36 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Sat 14 Feb 2009 - 18:30:23 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-devel. Please read the posting guide before posting to the list.

list of date sections of archive