zdata <- data.frame(x = runif(nn <- 2000))
zdata <- transform(zdata, p0 = logit(-1 + 2*x, inverse=TRUE),
y1 = rposnegbin(nn, munb=exp(0+2*x), size=exp(1)),
y2 = rposnegbin(nn, munb=exp(1+2*x), size=exp(1)))
zdata <- transform(zdata, y1 = ifelse(runif(nn) < p0, 0, y1),
y2 = ifelse(runif(nn) < p0, 0, y2))
with(zdata, table(y1))
with(zdata, table(y2))
fit <- vglm(cbind(y1,y2) ~ x, zanegbinomial, zdata, trace=TRUE)
coef(fit, matrix=TRUE)
head(fitted(fit))
head(predict(fit))
Run the code above in your browser using DataLab