x = runif(n <- 1000)
phi = logit(-0.5+1*x, inverse=TRUE)
munb = exp(3+x)
k = exp(0+2*x)
y1 = rzinb(n, phi, mu=munb, size=k)
y2 = rzinb(n, phi, mu=munb, size=k)
table(y1)["0"] / sum(table(y1))
table(y2)["0"] / sum(table(y2))
fit = vglm(cbind(y1,y2) ~ x, zinegbinomial(zero=NULL), trace=TRUE)
coef(fit, matrix=TRUE)
summary(fit)
cbind(fitted(fit), (1-phi) * munb)[1:5,]
vcov(fit)
Run the code above in your browser using DataLab