[R] Fit Negbin glm model with autoregressive correlation structure
Hello, I am attempting to estimate the effect of various variables on the time-series of counts of reported cattle stillbirths. We investigate the effect of day-of-week, month, holidays etc...and also the effect of non-temporal variables. We performed model comparisons between Gaussian glm, Poisson glm and negbin glm and the latter seems most appropriate for our data. We found that the residuals from our best model are not i.i.d. but follow an autoregressive process of order 5 , AR (5). I therefore wish to re-run this model after adding an AR(5) correlation structure in order to get unbiased estimates and standard errors for the variables retained in the model. In the past, I have been faced with a similar situation for a Gaussian glm and used the gls function with a corStruct object describing the within-group correlation structure. However, this would not work with our negbin model. Looking around on various help forums, I came across the possibilities of using generalized estimating equations instead. The gee function (in gee package) has a corstr object which would allow me to specify an AR process of whichever order but there is no option to include a negbin family. The geeglm function (in geepack package) does recognized the negbin family but only gives an option to fit an AR(1) correlation structure. The corstr object seems to have a userdefined option but it is unclear how it could be defined for an AR(5) process. In short, my questions are: *is it possible to include an AR (p) correlation structure directly into a negbin glm? How? * if GEE are the way forward, how can the the corstr object in geeglm be defined for an AR(p) process? Thank you for your suggestions, Dr Flavie Vial Veterinary Public Health Institute DCR-VPH, Vetsuisse Fakultät Schwarzenburgstrasse 155 CH-3003 Bern Switzerland flavie.v...@vetsuisse.unibe.chmailto:flavie.v...@vetsuisse.unibe.ch [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[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]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
