mydat = data.frame(x = seq(0, 1, len=(nsites <- 2000)))
mydat = transform(mydat, eta1 = -3 + 5 * x,
phi1 = logit(-1, inverse=TRUE),
oratio = exp(2))
mydat = transform(mydat, mu12 = cloglog(eta1, inverse=TRUE) * (1-phi1))
tmat = with(mydat, rbinom2.or(nsites, mu1=mu12, oratio=oratio, exch=TRUE))
mydat = transform(mydat, ybin1 = tmat[,1], ybin2 = tmat[,2])
with(mydat, table(ybin1,ybin2)) / nsites # For interest only
# Various plots of the data, for interest only
par(mfrow=c(2,2))
plot(jitter(ybin1) ~ x, data = mydat, col="blue")
plot(jitter(ybin2) ~ jitter(ybin1), data = mydat, col="blue")
plot(mu12 ~ x, data = mydat, col="blue", type="l", ylim=0:1,
ylab="Probability", main="Marginal probability and phi")
with(mydat, abline(h=phi1[1], col="red", lty="dashed"))
tmat2 = with(mydat, dbinom2.or(mu1=mu12, oratio=oratio, exch=TRUE))
with(mydat, matplot(x, tmat2, col=1:4, type="l", ylim=0:1,
ylab="Probability", main="Joint probabilities"))
# Now fit the model to the data.
fit = vglm(cbind(ybin1,ybin2) ~ x, fam=zipebcom, dat=mydat, trace=TRUE)
coef(fit, matrix=TRUE)
summary(fit)
vcov(fit)
Run the code above in your browser using DataLab