\HeaderA{simulate}{Simulate Responses}{simulate}
\keyword{models}{simulate}
\keyword{datagen}{simulate}
\begin{Description}\relax
Simulate one or more responses from the distribution
corresponding to a fitted model object.
\end{Description}
\begin{Usage}
\begin{verbatim}
simulate(object, nsim, seed, ...)
\end{verbatim}
\end{Usage}
\begin{Arguments}
\begin{ldescription}
\item[\code{object}] an object representing a fitted model.
\item[\code{nsim}] number of response vectors to simulate.  Defaults to \code{1}.
\item[\code{seed}] an object specifying if and how the random number
generator should be initialized (\sQuote{seeded}).\\
For the "lm" method, either \code{NULL} or an integer that will be
used in a call to \code{set.seed} before simulating the response
vectors.  If set, the value is saved as the \code{"seed"} attribute
of the returned value.  The default, \code{NULL} will not change the
random generator state, and return \code{\LinkA{.Random.seed}{.Random.seed}} as the
\code{"seed"} attribute, see \sQuote{Value}.

\item[\code{...}] additional optional arguments.
\end{ldescription}
\end{Arguments}
\begin{Details}\relax
This is a generic function.  Consult the individual modeling functions
for details on how to use this function.

Package \pkg{stats} has a method for \code{"lm"} objects which is used
for \code{\LinkA{lm}{lm}} and \code{\LinkA{glm}{glm}} fits.  There is a method
for fits from \code{glm.nb} in package \pkg{MASS}, and hence the case
of negative binomial families is not covered by the \code{"lm"}
method.

The methods for linear models fitted by \code{lm} or \code{glm(family
  = "gaussian")} assume that any weights which have been supplied are
inversely proportional to the error variance.  For other GLMs the
(optional) \code{simulate} component of the \code{\LinkA{family}{family}}
object is used---there is no appropriate simulation method for
\sQuote{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 \eqn{w_i}{} can be interpreted as meaning that
observation \eqn{i}{} is an average of \eqn{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 \code{\LinkA{gamma.shape}{gamma.shape}} in package
\pkg{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 \eqn{IG(\mu_i,
\lambda w_i)}{} (see
\url{http://en.wikipedia.org/wiki/Inverse_Gaussian_distribution})
where \eqn{lambda}{} is estimated by the inverse of the dispersion
estimate for the fit.  The variance is \eqn{\mu_i^3/(\lambda w_i)}{} and
hence inversely proportional to the prior weights.  The simulation is
done by function \code{\LinkA{rinvGaussian}{rinvGaussian}} from the
\pkg{SuppDists} package, which must be installed.
\end{Details}
\begin{Value}
Typically, a list of length \code{nsim} of simulated responses.  Where
appropriate the result can be a data frame (which is a special type of
list).

For the \code{"lm"} method, the result is a data frame with an
attribute \code{"seed"} containing the \code{seed} argument if not
\code{NULL} with \code{"kind"} attributes the value of
\code{as.list(\LinkA{RNGkind}{RNGkind}())}, otherwise (the default) the value of
\code{\LinkA{.Random.seed}{.Random.seed}} before the simulation was started.
\end{Value}
\begin{SeeAlso}\relax
\code{\LinkA{fitted.values}{fitted.values}} and \code{\LinkA{residuals}{residuals}} for related methods;
\code{\LinkA{glm}{glm}}, \code{\LinkA{lm}{lm}} for model fitting.

There are further examples in the \file{simulate.R} tests file in the
sources for package \pkg{stats}.
\end{SeeAlso}
\begin{Examples}
\begin{ExampleCode}
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
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
\end{ExampleCode}
\end{Examples}

