simulate | R Documentation |
Simulate one or more responses from the distribution corresponding to a fitted model object.
simulate(object, nsim = 1, seed = NULL, ...)
object |
an object representing a fitted model. |
nsim |
number of response vectors to simulate. Defaults to |
seed |
an object specifying if and how the random number
generator should be initialized (‘seeded’). |
... |
additional optional arguments. |
This is a generic function. Consult the individual modeling functions for details on how to use this function.
Package stats has a method for "lm"
objects which is used
for lm
and glm
fits. There is a method
for fits from glm.nb
in package MASS, and hence the
case of negative binomial families is not covered by the "lm"
method.
The methods for linear models fitted by lm
or glm(family
= "gaussian")
assume that any weights which have been supplied are
inversely proportional to the error variance. For other GLMs the
(optional) simulate
component of the family
object is used—there is no appropriate simulation method for
‘quasi’ models as they are specified only up to two moments.
For binomial and Poisson GLMs the dispersion is fixed at one. Integer prior weights w_i can be interpreted as meaning that observation i is an average of w_i observations, which is natural for binomials specified as proportions but less so for a Poisson, for which prior weights are ignored with a warning.
For a gamma GLM the shape parameter is estimated by maximum likelihood
(using function gamma.shape
in package
MASS). The interpretation of weights is as multipliers to a
basic shape parameter, since dispersion is inversely proportional to
shape.
For an inverse gaussian GLM the model assumed is
IG(μ_i, λ w_i) (see
https://en.wikipedia.org/wiki/Inverse_Gaussian_distribution)
where λ is estimated by the inverse of the dispersion
estimate for the fit. The variance is
μ_i^3/(λ w_i) and
hence inversely proportional to the prior weights. The simulation is
done by function rinvGauss
from the
SuppDists package, which must be installed.
Typically, a list of length nsim
of simulated responses. Where
appropriate the result can be a data frame (which is a special type of
list).
For the "lm"
method, the result is a data frame with an
attribute "seed"
. If argument seed
is NULL
, the
attribute is the value of .Random.seed
before the
simulation was started; otherwise it is the value of the argument with
a "kind"
attribute with value as.list(RNGkind())
.
RNG
about random number generation in R,
fitted.values
and residuals
for related methods;
glm
, lm
for model fitting.
There are further examples in the ‘simulate.R’ tests file in the sources for package stats.
x <- 1:5 mod1 <- lm(c(1:3, 7, 6) ~ x) S1 <- simulate(mod1, 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) # should be about the same as fitted(mod1) ## repeat identically: (sseed <- attr(S2, "seed")) # seed; RNGkind as attribute stopifnot(identical(S2, simulate(mod1, nsim = 200, seed = sseed))) ## To be sure about the proper RNGkind, e.g., after RNGversion("2.7.0") ## first set the RNG kind, then simulate do.call(RNGkind, attr(sseed, "kind")) identical(S2, simulate(mod1, nsim = 200, seed = sseed)) ## Binomial GLM examples yb1 <- matrix(c(4, 4, 5, 7, 8, 6, 6, 5, 3, 2), ncol = 2) modb1 <- glm(yb1 ~ x, family = binomial) S3 <- simulate(modb1, nsim = 4) # each column of S3 is a two-column matrix. x2 <- sort(runif(100)) yb2 <- rbinom(100, prob = plogis(2*(x2-1)), size = 1) yb2 <- factor(1 yb2, labels = c("failure", "success")) modb2 <- glm(yb2 ~ x2, family = binomial) S4 <- simulate(modb2, nsim = 4) # each column of S4 is a factor
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.