# Example 1. Dobson (1990) Page 93: Randomized Controlled Trial :
counts = c(18,17,15,20,10,20,25,13,12)
outcome = gl(3,1,9)
treatment = gl(3,3)
print(d.AD <- data.frame(treatment, outcome, counts))
vglm.D93 = vglm(counts ~ outcome + treatment, family=poissonff)
summary(vglm.D93)
# Example 2. Multinomial logit model
data(pneumo)
pneumo = transform(pneumo, let=log(exposure.time))
vglm(cbind(normal, mild, severe) ~ let, multinomial, pneumo)
# Example 3. Proportional odds model
fit = vglm(cbind(normal,mild,severe) ~ let, cumulative(par=TRUE), pneumo)
coef(fit, matrix=TRUE)
constraints(fit)
fit@x # LM model matrix
model.matrix(fit) # Larger VGLM model matrix
# Example 4. Bivariate logistic model
data(coalminers)
fit = vglm(cbind(nBnW, nBW, BnW, BW) ~ age, binom2.or, coalminers, trace=TRUE)
coef(fit, matrix=TRUE)
fit@y
# Example 5. The use of the xij argument
n = 1000
eyes = data.frame(lop = runif(n), rop = runif(n))
eyes = transform(eyes,
leye = ifelse(runif(n) < logit(-1+2*lop, inverse=TRUE), 1, 0),
reye = ifelse(runif(n) < logit(-1+2*rop, inverse=TRUE), 1, 0))
fit = vglm(cbind(leye,reye) ~ lop + rop + fill(lop),
binom2.or(exchangeable=TRUE, zero=3),
xij = op ~ lop + rop + fill(lop), data=eyes)
coef(fit)
coef(fit, matrix=TRUE)
coef(fit, matrix=TRUE, compress=FALSE)
# Here's one method to handle the xij argument with a term that
# produces more than one column in the model matrix.
POLY3 = function(x, ...) {
# A cubic
poly(c(x,...), 3)[1:length(x),]
}
fit = vglm(cbind(leye,reye) ~ POLY3(lop,rop) + POLY3(rop,lop) + fill(POLY3(lop,rop)),
binom2.or(exchangeable=TRUE, zero=3), data=eyes,
xij = POLY3(op) ~ POLY3(lop,rop) + POLY3(rop,lop) +
fill(POLY3(lop,rop)))
coef(fit)
coef(fit, matrix=TRUE)
coef(fit, matrix=TRUE, compress=FALSE)
predict(fit)[1:4,]
Run the code above in your browser using DataLab