simulate.vlm: Simulate Responses for VGLMs and VGAMs

In VGAM: Vector Generalized Linear and Additive Models

Simulate Responses for VGLMs and VGAMs

Description

Simulate one or more responses from the distribution corresponding to a fitted model object.

Usage

## S3 method for class 'vlm'
simulate(object, nsim = 1, seed = NULL, ...)

Arguments

object	an object representing a fitted model. Usually an object of class vglm-class or vgam-class.
nsim, seed	Same as simulate.
...	additional optional arguments.

Details

This is a methods function for simulate and hopefully should behave in a very similar manner. Only VGAM family functions with a simslot slot have been implemented for simulate.

Value

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.

Warning

With multiple response and/or multivariate responses, the order of the elements may differ. For some VGAM families, the order is n \times N \times F, where n is the sample size, N is nsim and F is ncol(fitted(vglmObject)). For other VGAM families, the order is n \times F \times N. An example of each is given below.

See Also

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. Also, categorical family functions: acat, cratio, sratio, cumulative, multinomial.

See also RNG about random number generation in R, vglm, vgam for model fitting.

Examples

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 0s

VGAM documentation built on Sept. 18, 2024, 9:09 a.m.