data(coalminers)
coalminers = transform(coalminers, age=(age-42)/5)
# Get the n x 4 matrix of counts
temp = vglm(cbind(nBnW,nBW,BnW,BW) ~ age, binom2.or, coalminers)
counts = round(c(weights(temp, type="prior")) * temp@y)
# Create a n x 2 matrix response for loglinb2()
fred = matrix(c(0,0, 0,1, 1,0, 1,1), 4, 2, byrow=TRUE)
yy = kronecker(matrix(1, nrow(counts), 1), fred)
wt = c(t(counts))
age = rep(coalminers$age, rep(4, length(coalminers$age)))
yy = yy[wt>0,]
age = age[wt>0]
wt = wt[wt>0]
fit = vglm(yy ~ age, loglinb2, trace=TRUE, wei=wt)
coef(fit, mat=TRUE) # Same! (at least for the log odds-ratio)
summary(fit)
# Try reconcile this with McCullagh and Nelder (1989), p.234
(0.166-0.131) / 0.027458 # 1.275 is approximately 1.25
Run the code above in your browser using DataLab