flavie.vial at vetsuisse.unibe.ch
2013-Aug-20 17:53 UTC
[R] Generating random values from zero inflated negative binomial distribution with rZINBI {gamlss.dist}
Good afternoon,
I have two years of count data (daily time-series data). My Y variable can be
modeled using a ZINB model with "day of the week" and
"month" effects (both x variables retained in the logistic and
negative binomial part of the model). Now I would like to simulate 100 years of
count data with the same characteristics as the 2 years of observations I have,
i.e. I want the predicted value for each day of the year to be set as the mean
of the ZINB distribution and to sample that distribution randomly to determine
the value for that day for a given year for each of 100 simulated years.
The function rZINBI from the {gamlss.dist} package sounds like what I need:
rZINBI(n, mu, sigma, nu)
I need to provide 3 parameters to the function to draw n random values from the
distribution.
mu
vector of positive means
sigma
vector of positive despersion parameter
nu
vector of zero probability parameter
Now, my question is how can I derive mu, sigma and nu from the output of my ZINB
model? Is mu the predicted mean count for a particular day-of-week*month
combination? Is sigma the overdispersion parameter (theta) from my ZINB model
(it appears from model output from example below that it is not)? What is nu?
######################################################
To give a reproducible example, I generate some fake data using rZINBI
#create data
count<- rZINBI(480, mu=300, sigma=10,nu=0.1)
day<-rep(1:5,96)# 5 days in week as week-ends are not included
month<-rep(rep(1:12,each=20),2)#let's assume 20 days in a month
data<-data.frame(count,day,month)
attach(data)
day<-factor(day)
month<-factor(month)
#activate library
library(pscl)
#fit ZINB model to data
f1d<-formula(count~day+month|day+month)
Nb1d<-zeroinfl(f1d,dist="negbin",link="logit")
summary(Nb1d)
######################################################
How can I extract the relevant information from my ZINB model (Nb1d) to feed
into the rZINBI function (knowing that I will have to do 12*5 times*100 for each
day*month combination for each 100 simulated years)?
Thank you for your suggestions
Dr Flavie Vial
Veterinary Public Health Institute
DCR-VPH, Vetsuisse Fakultät
Schwarzenburgstrasse 155
CH-3003 Bern
Switzerland
[[alternative HTML version deleted]]