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 bol...@ufl.edu / www.zoology.ufl.edu/bolker GPG key: www.zoology.ufl.edu/bolker/benbolker-publickey.asc
signature.asc
Description: OpenPGP digital signature
______________________________________________ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel