proposed simulate.glm method
Yes, thanks to Ben for getting the ball rolling. His code was more
streamlined than mine, pointing to further simplifications which I've
included in the extended version below.
The code for the additional families uses functions from MASS and
SuppDists - I wasn't sure about the best way to do this, so have just
used :: for now.
It appears to be working happily for both glm and gnm objects (no
gnm-specific code used).
Best wishes,
Heather
simulate.glm <- function (object, nsim = 1, seed = NULL, ...)
{
fam <- object$family$family
if(fam == "gaussian")
return(simulate.lm(object, nsim=nsim, seed=seed, ...))
if(!exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE))
runif(1) # initialize the RNG if necessary
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 <- object$fitted
ntot <- length(pred) * nsim
val <- switch(fam,
"poisson" = rpois(ntot, pred),
"binomial" = {
wts <- object$prior.weights
if (any(wts %% 1 != 0))
stop("cannot simulate from non-integer
prior.weights")
rbinom(ntot, size = wts, prob = pred)/wts
},
"Gamma" = {
shape <- MASS::gamma.shape(object)$alpha
rgamma(ntot, shape = shape, rate = shape/pred)
},
"inverse.gaussian" = {
lambda <- 1/summary(object)$dispersion
SuppDists::rinvGauss(ntot, nu = pred, lambda = lambda)
},
stop("family '", fam, "' not yet implemented"))
ans <- as.data.frame(matrix(val, ncol = nsim))
attr(ans, "seed") <- RNGstate
ans
}
Martin Maechler wrote:
Thanks a lot, Heather,
"HT" == Heather Turner <Heather.Turner at warwick.ac.uk>
on Fri, 13 Feb 2009 11:49:06 +0000 writes:
HT> Dear Martin,
HT> I think a simulate.glm method ought to be able to work for gnm objects
HT> too. David Firth and I started to work on this a long time ago, but
HT> stopped part-way through when simulate.lm was introduced, thinking that
HT> simulate.glm was probably in the pipeline and we were duplicating
HT> effort. Obviously we have let this slip when a contribution might have
HT> been useful. We developed a prototype for poisson, binomial, gaussian,
HT> gamma and inverse gaussian models which might be usefully merged with
HT> Ben's proposed simulate.glm. What's the best way to go about this? I
HT> would also like to test the proposed simulate.glm to check whether it
HT> will work with gnm objects or whether a simulate.gnm will be necessary.
In the mean time, private e-mail communications have started on
the subject, and yes, we are very insterested in finding
``the best'' possible way, probably making use of
Heather+David's code together with Ben's.
One alternative (not mentioned yet on R-devel), we've been
considering is to use simulate.lm() to also deal with "glm" (and
possibly "gnm") objects ``in one place''.
Martin
HT> Martin Maechler wrote:
>>>>>>> "BB" == Ben Bolker <bolker at ufl.edu>
>>>>>>> on Thu, 12 Feb 2009 11:29:14 -0500 writes:
>>
BB> I have found the "simulate" method (incorporated
BB> in some packages) very handy. As far as I can tell the
BB> only class for which simulate is actually implemented
BB> in base R is lm ... this is actually a little dangerous
BB> for a naive user who might be tempted to try
BB> simulate(X) where X is a glm fit instead, because
BB> it defaults to simulate.lm (since glm inherits from
BB> the lm class), and the answers make no sense ...
>>
BB> Here is my simulate.glm(), which is modeled on
BB> simulate.lm . It implements simulation for poisson
BB> and binomial (binary or non-binary) models, should
BB> be easy to implement others if that seems necessary.
>>
BB> I hereby request comments and suggest that it wouldn't
BB> hurt to incorporate it into base R ... (I will write
BB> docs for it if necessary, perhaps by modifying ?simulate --
BB> there is no specific documentation for simulate.lm)
>>
BB> cheers
BB> Ben Bolker
>>
>> [...............]
>>
>> Hi Ben,
>> thank you for your proposals.
>>
>> I agree that simulate.glm() has been in missing in some way,
>> till now, in particular, as, as you note, "glm" objects extend
>> "lm" ones and hence simulate(<glm>, ...) currently dispatches to
>> calling simulate.lm(....) which is only correct in the case of
>> the gaussian family.
>>
>> I have looked at your proposal a bit, already "improved" the
>> code slightly (e.g. re-include the comment you lost when you
>> ``copied'' simulate.lm(): In such cases, please work from the
>> source, not from what you get by print()ing
>> stats:::simulate.lm --- the source is either a recent tarball,
>> or the SVN repository, in this case, file
>> https://svn.r-project.org/R/trunk/src/library/stats/R/lm.R ]
>> and am planning to look at your and some own examples;
>> all with the goal to indeed include this in the R standard
>> 'stats' package in R-devel [to become R 2.9.0 in the future].
>>
>> About the help page: At the moment, I think that only a few
>> words would need to be added to the simulate help page,
>> i.e., https://svn.r-project.org/R/trunk/src/library/stats/man/simulate.Rd
>> and will be happy to receive a patch against this file.
>>
>> Thank you again, and best regards,
>> Martin Maechler, ETH Zurich