Re: [Rd] proposed simulate.glm method
On Sat, 14 Feb 2009, Nicholas Lewin-Koh wrote:
Well, my question wasn't that clear :-), but yes you mostly answered
it. I guess the one case I would be concerned is in Heather's code,
where the distribution to simulate from is chosen, that seemed to be
hard coded.
Rather, all known glm families in R that correspond to actual
probability distributions are listed.
So if I built a family object, say for a model that assumes errors
from a zipf distribution,
Hmm, plese explain how you get that into the GLM framework -- it is
pretty restrictive.
and I did have a predict method (which is a fair assumption) would that
fail because the rzipf function would not be accessed?
glm has a predict method, so why do you need one? Families do not
create additional classes of fit objects.
We could extend the definition of a family to have a 'simulate'
element, but then existing user-contributed families (principally the
negative binomial) would not have one and so this would not solve the
problem .
If you know of an actual R implementation of another glm family that
looks generally useful we'll consider adding it (but it seems that the
end user could also do so rather easily).
Nicholas
On Sat, 14 Feb 2009 20:10:26 +0100, "Martin Maechler"
said:
"NicLK" == Nicholas Lewin-Koh
on Sat, 14 Feb 2009 08:34:45 -0800 writes:
NicLK> Hi, For extended glms such as gams, gnm or other
NicLK> distributions such as negative binomial, would there
NicLK> need to be a separate simulate method?
Not necessarily, as I said, the "glm"s are now also dealt with
in simulate.lm() and Heather more or less confirmed that this
gives correct results for "gnm" objects.
For gam(), I'd strongly expect the same to apply, but there
maybe sophisticated gam() models where the result is currently
not correct. That's, BTW, also true for
simulate(lm(, weights), ...)
NicLK> Or, could the current framework, rather than
NicLK> stopping with an error look for the appropriate model
NicLK> matrix, coefficients, distribution function and
NicLK> family object to simulate from?
What do you mean?
A situation where there's no supported 'family'
or a situation where predict() does not work as it's
supposed in the current framework,
or
If there are such cases, we'd have to consider them together
with the corresponding package author. It may often make sense
fthen that the author changes his methods {predict(), ..} such
that the (now) extended simulate.lm() will work automatically;
Alternatively, the author can provide simulate.().
But I'm not sure I'm answering the question you've asked..
Martin
NicLK> Nicholas
>> Message: 9 Date: Fri, 13 Feb 2009 21:27:57 +0100 From:
>> Martin Maechler Subject: Re:
>> [Rd] proposed simulate.glm method To: Heather Turner
>> Cc: r-devel@r-project.org,
>> Martin Maechler Message-ID:
>> <18837.55245.15158.29...@cmath-5.math.ethz.ch>
>> Content-Type: text/plain; charset=us-ascii
>>
>> Thank you, Heather and Ben,
>>
>>>>>>> "HT" == Heather Turner
>> >>>>> on Fri, 13 Feb 2009
>> 15:52:37 + writes:
>>
HT> Yes, thanks to Ben for getting the ball rolling. His
HT> code was more streamlined than mine, pointing to further
HT> simplifications which I've included in the extended
HT> version below.
>>
HT> The code for the additional families uses functions from
HT> MASS and SuppDists - I wasn't sure about the best way to
HT> do this, so have just used :: for now.
>>
HT> It appears to be working happily for both glm and gnm
HT> objects (no gnm-specific code used).
>>
HT> Best wishes,
>>
HT> Heather
>>
>> []
>>
>> I have now followed Brian Ripley's suggetion to just
>> extend simulate.lm() to also deal with "glm" objects, but
>> using Heather's suggestions for the different families;
>> I've just commited src/library/stats/R/lm.R with the new
>> code. (get it from svn.r-project.org/R/trunk/ or this
>> night's R-devel tarball).
>>
>> One difference to your propsal: Instead of just
>> object$fitted , the code is using fitted(object)
>> ... something which should properly to the na.action
>> used.
>>
>> For the (MASS and) SuppDists package requirement, I'm
>> using the following
>>
>> if(is.null(tryCatch(loadNamespace("SuppDists"), error =
>> function(e) NULL))) stop("Need CRAN package '
Re: [Rd] proposed simulate.glm method
Hi,
Well, my question wasn't that clear :-), but yes you mostly answered it.
I guess
the one case I would be concerned is in Heather's code, where the
distribution
to simulate from is chosen, that seemed to be hard coded. So if
I built a family object, say for a model that assumes errors from a zipf
distribution,
and I did have a predict method (which is a fair assumption) would that
fail because the rzipf function would not be accessed?
Nicholas
On Sat, 14 Feb 2009 20:10:26 +0100, "Martin Maechler"
said:
> >>>>> "NicLK" == Nicholas Lewin-Koh
> >>>>> on Sat, 14 Feb 2009 08:34:45 -0800 writes:
>
> NicLK> Hi, For extended glms such as gams, gnm or other
> NicLK> distributions such as negative binomial, would there
> NicLK> need to be a separate simulate method?
>
> Not necessarily, as I said, the "glm"s are now also dealt with
> in simulate.lm() and Heather more or less confirmed that this
> gives correct results for "gnm" objects.
>
> For gam(), I'd strongly expect the same to apply, but there
> maybe sophisticated gam() models where the result is currently
> not correct. That's, BTW, also true for
> simulate(lm(, weights), ...)
>
> NicLK> Or, could the current framework, rather than
> NicLK> stopping with an error look for the appropriate model
> NicLK> matrix, coefficients, distribution function and
> NicLK> family object to simulate from?
>
> What do you mean?
> A situation where there's no supported 'family'
> or a situation where predict() does not work as it's
> supposed in the current framework,
> or
>
> If there are such cases, we'd have to consider them together
> with the corresponding package author. It may often make sense
> fthen that the author changes his methods {predict(), ..} such
> that the (now) extended simulate.lm() will work automatically;
> Alternatively, the author can provide simulate.().
>
> But I'm not sure I'm answering the question you've asked..
> Martin
>
> NicLK> Nicholas
>
>
> >> Message: 9 Date: Fri, 13 Feb 2009 21:27:57 +0100 From:
> >> Martin Maechler Subject: Re:
> >> [Rd] proposed simulate.glm method To: Heather Turner
> >> Cc: r-devel@r-project.org,
> >> Martin Maechler Message-ID:
> >> <18837.55245.15158.29...@cmath-5.math.ethz.ch>
> >> Content-Type: text/plain; charset=us-ascii
> >>
> >> Thank you, Heather and Ben,
> >>
> >> >>>>> "HT" == Heather Turner
> >> >>>>> on Fri, 13 Feb 2009
> >> 15:52:37 + writes:
> >>
> HT> Yes, thanks to Ben for getting the ball rolling. His
> HT> code was more streamlined than mine, pointing to further
> HT> simplifications which I've included in the extended
> HT> version below.
> >>
> HT> The code for the additional families uses functions from
> HT> MASS and SuppDists - I wasn't sure about the best way to
> HT> do this, so have just used :: for now.
> >>
> HT> It appears to be working happily for both glm and gnm
> HT> objects (no gnm-specific code used).
> >>
> HT> Best wishes,
> >>
> HT> Heather
> >>
> >> []
> >>
> >> I have now followed Brian Ripley's suggetion to just
> >> extend simulate.lm() to also deal with "glm" objects, but
> >> using Heather's suggestions for the different families;
> >> I've just commited src/library/stats/R/lm.R with the new
> >> code. (get it from svn.r-project.org/R/trunk/ or this
> >> night's R-devel tarball).
> >>
> >> One difference to your propsal: Instead of just
> >> object$fitted , the code is using fitted(object)
> >> ... something which should properly to the na.action
> >> used.
> >>
> >> For the (MASS and) SuppDists package requirement, I'm
> >> using the following
> >>
> >> if(is.null(tryCatch(loadNamespace("SuppDists"), error =
> >> function(e) NULL))) stop("Need CRAN package 'SuppDists'
> >> for 'inverse.gaussian' family")
> >>
> >>
> >> I've not yet updated the help page for simulate(), and
> >> have only tested relatively few cases for binomial,
> >> poisson and Gamma. I've wanted to expose this to you, so
> >> you can provide more feedback and possibly even a patch
> >> to
> >> svn.r-project.org/R/trunk/src/library/stats/man/simulate.Rd
> >>
> >> Martin
> >>
> >>
> >>
>
> __
> NicLK> R-devel@r-project.org mailing list
> NicLK> https://stat.ethz.ch/mailman/listinfo/r-devel
__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] proposed simulate.glm method
>>>>> "NicLK" == Nicholas Lewin-Koh
>>>>> on Sat, 14 Feb 2009 08:34:45 -0800 writes:
NicLK> Hi, For extended glms such as gams, gnm or other
NicLK> distributions such as negative binomial, would there
NicLK> need to be a separate simulate method?
Not necessarily, as I said, the "glm"s are now also dealt with
in simulate.lm() and Heather more or less confirmed that this
gives correct results for "gnm" objects.
For gam(), I'd strongly expect the same to apply, but there
maybe sophisticated gam() models where the result is currently
not correct. That's, BTW, also true for
simulate(lm(, weights), ...)
NicLK> Or, could the current framework, rather than
NicLK> stopping with an error look for the appropriate model
NicLK> matrix, coefficients, distribution function and
NicLK> family object to simulate from?
What do you mean?
A situation where there's no supported 'family'
or a situation where predict() does not work as it's
supposed in the current framework,
or
If there are such cases, we'd have to consider them together
with the corresponding package author. It may often make sense
fthen that the author changes his methods {predict(), ..} such
that the (now) extended simulate.lm() will work automatically;
Alternatively, the author can provide simulate.().
But I'm not sure I'm answering the question you've asked..
Martin
NicLK> Nicholas
>> Message: 9 Date: Fri, 13 Feb 2009 21:27:57 +0100 From:
>> Martin Maechler Subject: Re:
>> [Rd] proposed simulate.glm method To: Heather Turner
>> Cc: r-devel@r-project.org,
>> Martin Maechler Message-ID:
>> <18837.55245.15158.29...@cmath-5.math.ethz.ch>
>> Content-Type: text/plain; charset=us-ascii
>>
>> Thank you, Heather and Ben,
>>
>> >>>>> "HT" == Heather Turner
>> >>>>> on Fri, 13 Feb 2009
>> 15:52:37 + writes:
>>
HT> Yes, thanks to Ben for getting the ball rolling. His
HT> code was more streamlined than mine, pointing to further
HT> simplifications which I've included in the extended
HT> version below.
>>
HT> The code for the additional families uses functions from
HT> MASS and SuppDists - I wasn't sure about the best way to
HT> do this, so have just used :: for now.
>>
HT> It appears to be working happily for both glm and gnm
HT> objects (no gnm-specific code used).
>>
HT> Best wishes,
>>
HT> Heather
>>
>> []
>>
>> I have now followed Brian Ripley's suggetion to just
>> extend simulate.lm() to also deal with "glm" objects, but
>> using Heather's suggestions for the different families;
>> I've just commited src/library/stats/R/lm.R with the new
>> code. (get it from svn.r-project.org/R/trunk/ or this
>> night's R-devel tarball).
>>
>> One difference to your propsal: Instead of just
>> object$fitted , the code is using fitted(object)
>> ... something which should properly to the na.action
>> used.
>>
>> For the (MASS and) SuppDists package requirement, I'm
>> using the following
>>
>> if(is.null(tryCatch(loadNamespace("SuppDists"), error =
>> function(e) NULL))) stop("Need CRAN package 'SuppDists'
>> for 'inverse.gaussian' family")
>>
>>
>> I've not yet updated the help page for simulate(), and
>> have only tested relatively few cases for binomial,
>> poisson and Gamma. I've wanted to expose this to you, so
>> you can provide more feedback and possibly even a patch
>> to
>> svn.r-project.org/R/trunk/src/library/stats/man/simulate.Rd
>>
>> Martin
>>
>>
>>
__
NicLK> R-devel@r-project.org mailing list
NicLK> https://stat.ethz.ch/mailman/listinfo/r-devel
__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] proposed simulate.glm method
Hi,
For extended glms such as gams, gnm or other distributions
such as negative binomial, would there need to be a separate simulate
method?
Or, could the current framework, rather than stopping with an error
look for the appropriate model matrix, coefficients, distribution
function and family object
to simulate from?
Nicholas
> Message: 9
> Date: Fri, 13 Feb 2009 21:27:57 +0100
> From: Martin Maechler
> Subject: Re: [Rd] proposed simulate.glm method
> To: Heather Turner
> Cc: r-devel@r-project.org, Martin Maechler
>
> Message-ID: <18837.55245.15158.29...@cmath-5.math.ethz.ch>
> Content-Type: text/plain; charset=us-ascii
>
> Thank you, Heather and Ben,
>
> >>>>> "HT" == Heather Turner
> >>>>> on Fri, 13 Feb 2009 15:52:37 + writes:
>
> HT> Yes, thanks to Ben for getting the ball rolling. His
> HT> code was more streamlined than mine, pointing to further
> HT> simplifications which I've included in the extended
> HT> version below.
>
> HT> The code for the additional families uses functions from
> HT> MASS and SuppDists - I wasn't sure about the best way to
> HT> do this, so have just used :: for now.
>
> HT> It appears to be working happily for both glm and gnm
> HT> objects (no gnm-specific code used).
>
> HT> Best wishes,
>
> HT> Heather
>
> []
>
> I have now followed Brian Ripley's suggetion to just extend
> simulate.lm() to also deal with "glm" objects, but using
> Heather's suggestions for the different families;
> I've just commited src/library/stats/R/lm.R with the new code.
> (get it from svn.r-project.org/R/trunk/ or this night's R-devel
> tarball).
>
> One difference to your propsal: Instead of just
> object$fitted , the code is using
> fitted(object) ... something which should properly to the na.action
> used.
>
> For the (MASS and) SuppDists package requirement, I'm using
> the following
>
> if(is.null(tryCatch(loadNamespace("SuppDists"),
> error = function(e) NULL)))
> stop("Need CRAN package 'SuppDists' for 'inverse.gaussian' family")
>
>
> I've not yet updated the help page for simulate(),
> and have only tested relatively few cases for binomial, poisson
> and Gamma.
> I've wanted to expose this to you, so you can provide more
> feedback and possibly even a patch to
>svn.r-project.org/R/trunk/src/library/stats/man/simulate.Rd
>
> Martin
>
>
>
__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] proposed simulate.glm method
Thank you, Heather and Ben,
> "HT" == Heather Turner
> on Fri, 13 Feb 2009 15:52:37 + writes:
HT> Yes, thanks to Ben for getting the ball rolling. His
HT> code was more streamlined than mine, pointing to further
HT> simplifications which I've included in the extended
HT> version below.
HT> The code for the additional families uses functions from
HT> MASS and SuppDists - I wasn't sure about the best way to
HT> do this, so have just used :: for now.
HT> It appears to be working happily for both glm and gnm
HT> objects (no gnm-specific code used).
HT> Best wishes,
HT> Heather
[]
I have now followed Brian Ripley's suggetion to just extend
simulate.lm() to also deal with "glm" objects, but using
Heather's suggestions for the different families;
I've just commited src/library/stats/R/lm.R with the new code.
(get it from svn.r-project.org/R/trunk/ or this night's R-devel
tarball).
One difference to your propsal: Instead of just
object$fitted , the code is using
fitted(object) ... something which should properly to the na.action
used.
For the (MASS and) SuppDists package requirement, I'm using
the following
if(is.null(tryCatch(loadNamespace("SuppDists"),
error = function(e) NULL)))
stop("Need CRAN package 'SuppDists' for 'inverse.gaussian' family")
I've not yet updated the help page for simulate(),
and have only tested relatively few cases for binomial, poisson
and Gamma.
I've wanted to expose this to you, so you can provide more
feedback and possibly even a patch to
svn.r-project.org/R/trunk/src/library/stats/man/simulate.Rd
Martin
__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] proposed simulate.glm method
If you are generalizing this, the saving of the RNG information to reproduce normally distribution random number also needs to save the normal generator information. (This looks like an omission in simulate.lm.) You might want to consider adding the simple functions setRNG and getRNG from my setRNG package. Roughly, the functions make the first ten lines of simulate.lm or Ben's code into a function call, which may not be so important other than the missing normal information, but having a standardized object with all the information seems to be useful. The package also has documentation and tests, which can be helpful for other package builders that define simulate methods. It has been around for a long time so has been pretty well tested. Paul Gilbert Martin Maechler wrote: Thanks a lot, Heather, "HT" == Heather Turner on Fri, 13 Feb 2009 11:49:06 + 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 >>> 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(, ...) 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 __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel La version française suit le texte anglais.
Re: [Rd] 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
>> on Fri, 13 Feb 2009 11:49:06 + 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
> >>> 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(, ...) currently dispatches t
Re: [Rd] proposed simulate.glm method
Thanks a lot, Heather, > "HT" == Heather Turner > on Fri, 13 Feb 2009 11:49:06 + 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 >>> 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(, ...) 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 __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] proposed simulate.glm method
Dear Martin, I think a simulate.glm method ought to be able to work for gnm objects too. David Firth and I started to work on this a long time ago, but stopped part-way through when simulate.lm was introduced, thinking that simulate.glm was probably in the pipeline and we were duplicating effort. Obviously we have let this slip when a contribution might have been useful. We developed a prototype for poisson, binomial, gaussian, gamma and inverse gaussian models which might be usefully merged with Ben's proposed simulate.glm. What's the best way to go about this? I would also like to test the proposed simulate.glm to check whether it will work with gnm objects or whether a simulate.gnm will be necessary. Thanks and best regards, Heather Martin Maechler wrote: >> "BB" == Ben Bolker >> 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(, ...) 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 > > __ > R-devel@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-devel __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] proposed simulate.glm method
> "BB" == Ben Bolker > 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(, ...) 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 __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] proposed simulate.glm method
Elsewhere (at least in lme4), refit(sim(model)) does the same thing [and so one would need something like apply(sim(model,1000),2,refit)]. sim() is quite interesting, as is Zelig, but I'm not sure I am ready to leap to it yet -- this was basically a suggestion that simulate.glm could be included in "vanilla" R ... Also (for better or worse), it looks like sim() also does parametric bootstrapping on the parameter values, whereas simulate.[g]lm() just uses "plug-in" estimates. cheers Ben Bolker Alex D'Amour wrote: > There is functionality similar to this included in the Zelig package > with it's "sim" method. The "sim" method goes a step further and > replicates the fitted model's analysis on the generated datasets as > well. I would suggest taking a look -- Zelig supports most (if not > all) glm models and a wide range of others. > > The Zelig maintainers' site can be found at: http://gking.harvard.edu/zelig/. > > Full disclosure: I am an employee of the Institute for Quantitative > Social Science, which performs most of the development and support for > the Zelig package. > > Best, > Alex D'Amour > Statistical Programmer > Harvard Institute for Quantitative Social Science > > > 2009/2/12 Ben Bolker : >> 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 >> -- >> 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 >> >> >> __ >> R-devel@r-project.org mailing list >> https://stat.ethz.ch/mailman/listinfo/r-devel >> >> -- 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 __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] proposed simulate.glm method
There is functionality similar to this included in the Zelig package with it's "sim" method. The "sim" method goes a step further and replicates the fitted model's analysis on the generated datasets as well. I would suggest taking a look -- Zelig supports most (if not all) glm models and a wide range of others. The Zelig maintainers' site can be found at: http://gking.harvard.edu/zelig/. Full disclosure: I am an employee of the Institute for Quantitative Social Science, which performs most of the development and support for the Zelig package. Best, Alex D'Amour Statistical Programmer Harvard Institute for Quantitative Social Science 2009/2/12 Ben Bolker : > > 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 > -- > 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 > > > __ > R-devel@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-devel > > __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
