rnbinom.gf {spass}

Generate Time Series with Negative Binomial Distribution and Multivariate Gamma Frailty with Autoregressive Correlation Structure of Order One

Description

rnbinom.gf generates one or more independent time series following the Gamma frailty model. The generated data has negative binomial marginal distribution and the underlying multivariate Gamma frailty an autoregressive covariance structure.

Usage

rnbinom.gf(n, size, mu, rho, tp)

Arguments

n

number of observations. If length(n) > 1, the length is taken to be the number required.

size

dispersion parameter (the shape parameter of the gamma mixing distribution). Must be strictly positive, need not be integer.

mu

vector of means of time points: see 'Details'.

rho

correlation coefficient of the underlying autoregressive Gamma frailty. Must be between 0 and 1.

tp

number of observed time points.

Details

The generated marginal negative binomial distribution with mean mu = \mu and size = \eta has density

(\mu/(\mu+\eta))^x \Gamma(x + \eta)/(\Gamma(x+1)\Gamma(\eta)) (\eta/(\mu+\eta))^\eta

for 0 < \mu, 0 < \eta and x=0, 1, 2, .... Hereby, each entry of vector mu corresponds to one time point. Therefore, each timepoint can have its distinct mean.

Within the Gamma frailty model, the correlation between two frailties of time points t and s for rho = \rho is given by

\rho^|t-s|

for 0 \le \rho \le 1. Note: this does not correspond to the correlation of observations.

Value

rnbinom.gf returns a matrix of dimension n x tp with marginal negative binomial distribution with means mu, common dispersion parameter size and a correlation induce by the autoregressive multivariate Gamma frailty.

Source

rnbinom.gf computes observations from a Gamma frailty model by Fiocco et. al. 2009 using code contributed by Thomas Asendorf.

References

Fiocco M, Putter H, Van Houwelingen JC, (2009), A new serially correlated gamma-frailty process for longitudinal count data Biostatistics Vol. 10, No. 2, pp. 245-257.

Examples

set.seed(8)
random<-rnbinom.gf(n=1000, size=0.6, mu=1:6, rho=0.8, tp=6)
cor(random)

#Check the marginal distribution of time point 3
plot(table(random[,3])/1000, xlab="Probability", ylab="Observation")
lines(0:26, dnbinom(0:26, mu=3, size=0.6), col="red")
legend("topright",legend=c("Theoretical Marginal Distribution", "Observed Distribution"),
  col=c("red", "black"), lty=1, lwd=c(1,2))