spass (version 1.2)

rnbinom.gf: 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.

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.

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.

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

# NOT RUN {
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))

# }

Run the code above in your browser using DataCamp Workspace