[Rd] proposed simulate.glm method
Ben Bolker
bolker at ufl.edu
Thu Feb 12 17:29:14 CET 2009
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
-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 260 bytes
Desc: OpenPGP digital signature
URL: <https://stat.ethz.ch/pipermail/r-devel/attachments/20090212/fb93e658/attachment.bin>
More information about the R-devel
mailing list