# Fit the model in Table 6.7 in McCullagh and Nelder (1989)
data(coalminers)
coalminers = transform(coalminers, Age = (age - 42) / 5)
fit = vglm(cbind(nBnW,nBW,BnW,BW) ~ Age, binom2.or(zero=NULL), coalminers)
fitted(fit)
summary(fit)
coef(fit, matrix=TRUE)
with(coalminers, matplot(Age, fitted(fit), type="l", las=1,
xlab="(age - 42) / 5", lwd=2))
with(coalminers, matpoints(Age, fit@y, col=1:4))
legend(x=-4, y=0.5, lty=1:4, col=1:4, lwd=2,
legend=c("1 = (Breathlessness=0, Wheeze=0)",
"2 = (Breathlessness=0, Wheeze=1)",
"3 = (Breathlessness=1, Wheeze=0)",
"4 = (Breathlessness=1, Wheeze=1)"))
Run the code above in your browser using DataLab