nbdat = data.frame(x = runif(nn <- 1000))
nbdat = transform(nbdat, phi = logit(-0.5+1*x, inverse=TRUE),
munb = exp(3+x),
k = exp(0+2*x))
nbdat = transform(nbdat, y1 = rzinegbin(nn, phi, mu=munb, size=k),
y2 = rzinegbin(nn, phi, mu=munb, size=k))
with(nbdat, table(y1)["0"] / sum(table(y1)))
fit = vglm(cbind(y1,y2) ~ x, zinegbinomial(zero=NULL), nbdat, trace=TRUE)
coef(fit, matrix=TRUE)
summary(fit)
head(cbind(fitted(fit), with(nbdat, (1-phi) * munb)))
round(vcov(fit), 3)
Run the code above in your browser using DataLab