# NOT RUN {
# Fit the model in Table 6.7 in McCullagh and Nelder (1989)
coalminers <- transform(coalminers, Age = (age - 42) / 5)
fit <- vglm(cbind(nBnW, nBW, BnW, BW) ~ Age,
binom2.or(zero = NULL), data = coalminers)
fitted(fit)
summary(fit)
coef(fit, matrix = TRUE)
c(weights(fit, type = "prior")) * fitted(fit) # Table 6.8
# }
# NOT RUN {
with(coalminers, matplot(Age, fitted(fit), type = "l", las = 1,
xlab = "(age - 42) / 5", lwd = 2))
with(coalminers, matpoints(Age, depvar(fit), 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)"))
# }
# NOT RUN {
# Another model: pet ownership
# }
# NOT RUN {
data(xs.nz, package = "VGAMdata")
# More homogeneous:
petdata <- subset(xs.nz, ethnicity == "European" & age < 70 & sex == "M")
petdata <- na.omit(petdata[, c("cat", "dog", "age")])
summary(petdata)
with(petdata, table(cat, dog)) # Can compute the odds ratio
fit <- vgam(cbind((1-cat) * (1-dog), (1-cat) * dog,
cat * (1-dog), cat * dog) ~ s(age, df = 5),
binom2.or(zero = 3), data = petdata, trace = TRUE)
colSums(depvar(fit))
coef(fit, matrix = TRUE)
# }
# NOT RUN {
# }
# NOT RUN {
# Plot the estimated probabilities
ooo <- order(with(petdata, age))
matplot(with(petdata, age)[ooo], fitted(fit)[ooo, ], type = "l",
xlab = "Age", ylab = "Probability", main = "Pet ownership",
ylim = c(0, max(fitted(fit))), las = 1, lwd = 1.5)
legend("topleft", col=1:4, lty = 1:4, leg = c("no cat or dog ",
"dog only", "cat only", "cat and dog"), lwd = 1.5)
# }
Run the code above in your browser using DataLab