Simulate Responses for VGLMs and VGAMs
Simulate one or more responses from the distribution corresponding to a fitted model object.
## S3 method for class 'vlm' simulate(object, nsim = 1, seed = NULL, ...)
object |
an object representing a fitted model.
Usually an object of class
|
nsim, seed |
Same as |
... |
additional optional arguments. |
Similar to simulate.
Note that many VGAM family functions can handle multiple responses.
This can result in a longer data frame with more rows
(nsim multiplied by n rather than the
ordinary n).
In the future an argument may be available so that there
is always n rows no matter how many responses were
inputted.
With multiple response and/or multivariate responses,
the order of the elements may differ.
For some VGAM families, the order is
n x N x F,
where n is the sample size,
N is nsim and
F is ncol(fitted(vglmObject)).
For other VGAM families, the order is
n x F x N.
An example of each is given below.
Currently the VGAM family functions with a
simslot slot are:
alaplace1,
alaplace2,
betabinomial,
betabinomialff,
betaR,
betaff,
biamhcop,
bifrankcop,
bilogistic,
binomialff,
binormal,
binormalcop,
biclaytoncop,
cauchy,
cauchy1,
chisq,
dirichlet,
dagum,
erlang,
exponential,
bifgmcop,
fisk,
gamma1,
gamma2,
gammaR,
gengamma.stacy,
geometric,
gompertz,
gumbelII,
hzeta,
inv.lomax,
inv.paralogistic,
kumar,
lgamma1,
lgamma3,
lindley,
lino,
logff,
logistic1,
logistic,
lognormal,
lomax,
makeham,
negbinomial,
negbinomial.size,
paralogistic,
perks,
poissonff,
posnegbinomial,
posnormal,
pospoisson,
polya,
polyaR,
posbinomial,
rayleigh,
riceff,
simplex,
sinmad,
slash,
studentt,
studentt2,
studentt3,
triangle,
uninormal,
yulesimon,
zageometric,
zageometricff,
zanegbinomial,
zanegbinomialff,
zapoisson,
zapoissonff,
zigeometric,
zigeometricff,
zinegbinomial,
zipf,
zipoisson,
zipoissonff.
nn <- 10; mysize <- 20; set.seed(123)
bdata <- data.frame(x2 = rnorm(nn))
bdata <- transform(bdata,
y1 = rbinom(nn, size = mysize, p = logitlink(1+x2, inverse = TRUE)),
y2 = rbinom(nn, size = mysize, p = logitlink(1+x2, inverse = TRUE)),
f1 = factor(as.numeric(rbinom(nn, size = 1,
p = logitlink(1+x2, inverse = TRUE)))))
(fit1 <- vglm(cbind(y1, aaa = mysize - y1) ~ x2, # Matrix response (2-colns)
binomialff, data = bdata))
(fit2 <- vglm(f1 ~ x2, binomialff, model = TRUE, data = bdata)) # Factor response
set.seed(123); simulate(fit1, nsim = 8)
set.seed(123); c(simulate(fit2, nsim = 3)) # Use c() when model = TRUE
# An n x N x F example
set.seed(123); n <- 100
bdata <- data.frame(x2 = runif(n), x3 = runif(n))
bdata <- transform(bdata, y1 = rnorm(n, 1 + 2 * x2),
y2 = rnorm(n, 3 + 4 * x2))
fit1 <- vglm(cbind(y1, y2) ~ x2, binormal(eq.sd = TRUE), data = bdata)
nsim <- 1000 # Number of simulations for each observation
my.sims <- simulate(fit1, nsim = nsim)
dim(my.sims) # A data frame
aaa <- array(unlist(my.sims), c(n, nsim, ncol(fitted(fit1)))) # n by N by F
summary(rowMeans(aaa[, , 1]) - fitted(fit1)[, 1]) # Should be all 0s
summary(rowMeans(aaa[, , 2]) - fitted(fit1)[, 2]) # Should be all 0s
# An n x F x N example
n <- 100; set.seed(111); nsim <- 1000
zdata <- data.frame(x2 = runif(n))
zdata <- transform(zdata, lambda1 = loglink(-0.5 + 2 * x2, inverse = TRUE),
lambda2 = loglink( 0.5 + 2 * x2, inverse = TRUE),
pstr01 = logitlink( 0, inverse = TRUE),
pstr02 = logitlink(-1.0, inverse = TRUE))
zdata <- transform(zdata, y1 = rzipois(n, lambda = lambda1, pstr0 = pstr01),
y2 = rzipois(n, lambda = lambda2, pstr0 = pstr02))
zip.fit <- vglm(cbind(y1, y2) ~ x2, zipoissonff, data = zdata, crit = "coef")
my.sims <- simulate(zip.fit, nsim = nsim)
dim(my.sims) # A data frame
aaa <- array(unlist(my.sims), c(n, ncol(fitted(zip.fit)), nsim)) # n by F by N
summary(rowMeans(aaa[, 1, ]) - fitted(zip.fit)[, 1]) # Should be all 0s
summary(rowMeans(aaa[, 2, ]) - fitted(zip.fit)[, 2]) # Should be all 0sPlease choose more modern alternatives, such as Google Chrome or Mozilla Firefox.