size <- 10 # Number of trials; N in the notation above
nn <- 200
zibdata <- data.frame(pstr0 = logit( 0, inverse = TRUE), # 0.50
mubin = logit(-1, inverse = TRUE), # Mean of usual binomial
sv = rep(size, length = nn))
zibdata <- transform(zibdata,
y = rzibinom(nn, size = sv, prob = mubin, pstr0 = pstr0))
with(zibdata, table(y))
fit <- vglm(cbind(y, sv - y) ~ 1, zibinomial, zibdata, trace = TRUE)
fit <- vglm(cbind(y, sv - y) ~ 1, zibinomial, zibdata, trace = TRUE, stepsize = 0.5)
coef(fit, matrix = TRUE)
Coef(fit) # Useful for intercept-only models
fit@misc$pobs0 # Estimate of P(Y = 0)
head(fitted(fit))
with(zibdata, mean(y)) # Compare this with fitted(fit)
summary(fit)
Run the code above in your browser using DataLab