From: Ben Bolker <bolker_at_ufl.edu>

Date: Thu, 12 Feb 2009 11:29:14 -0500

ans

}

S1B <- simulate(mod2, nsim = 4)

## repeat the simulation:

.Random.seed <- attr(S1, "seed")

identical(S1, simulate(mod1, nsim = 4))

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.ascReceived on Thu 12 Feb 2009 - 15:50:36 GMT

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

- application/pgp-signature attachment: OpenPGP digital signature

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.
*