# Some R code to accompany
# "The VGAM Package for Categorical Data Analysis",
# by Yee, T. W. (2009)
# Journal of Statistical Software, 32.

# Last modified: 20091106


# ----------------------------------------------------------------------
# Load the package in

library("VGAM")
ps.options(pointsize = 12)
options(width = 72, digits = 4)


# ----------------------------------------------------------------------
# Bradley Terry code

journal <- c("Biometrika", "Comm.Statist", "JASA", "JRSS-B")
squaremat <- matrix(c( NA, 33, 320, 284,   730, NA, 813, 276,
             498, 68,  NA, 325,   221, 17, 142, NA), 4, 4)
dimnames(squaremat) <- list(winner = journal, loser = journal)
Brat(squaremat)


# ----------------------------------------------------------------------
# Genetic example

abodat <- data.frame(A = 725, B = 258, AB = 72, O = 1073)
fit <- vglm(cbind(A,B,AB,O) ~ 1, ABO, abodat)
coef(fit, matrix = TRUE)
Coef(fit) # Estimated pA and pB



# ----------------------------------------------------------------------
# 2008 World Fly Fishing Data example

head(wffc.nc, 5)

fnc <- transform(wffc.nc,
                finame = factor(iname),
                fsector = factor(sector),
                fday = factor(ceiling(session / 2)),
                mornaft = 1 - (session %% 2),
                fbeatboat = factor(beatboat))

fnc <- fnc[with(fnc, !is.element(comid, c(99,72,80,93,45,71,97,78))),] 
fnc <- transform(fnc,
                ordnum = ifelse(numbers <= 02, "few",
                         ifelse(numbers <= 10, "more", "most")))
fnc$ordnum <- ordered(fnc$ordnum, levels = c("few","more","most"))

with(fnc, table(ordnum))


# Fit a proportional odds model
fit.pom <-
  vglm(ordnum ~ fsector + mornaft + fday + finame,
       cumulative(parallel = TRUE, reverse = TRUE),
       data = fnc)



# Check that there are no problems with the input.
head(fit.pom@y, 3)        # Poor style
head(depvar(fit.pom), 3)  # Better style
colSums(depvar(fit.pom))

head(coef(fit.pom, matrix = TRUE), 10)


# Standard errors and Wald statistics may be obtained by
head(coef(summary(fit.pom)), 10)



# Check the proportional odds assumption with respect
# to the variable mornaft.
fit.ppom <-
  vglm(ordnum ~ fsector + mornaft + fday + finame,
       cumulative(parallel = FALSE ~ 1 + mornaft, reverse = TRUE),
       data = fnc)
head(coef(fit.ppom, mat = TRUE),  8)



# Likelihood ratio test $p$-value which is non-significant.
pchisq(deviance(fit.pom) - deviance(fit.ppom),
       df = df.residual(fit.pom) - df.residual(fit.ppom),
       lower.tail = FALSE)


# Final model
fit2.ppom <-
  vglm(ordnum ~ fsector + mornaft + fday + finame,
       cumulative(parallel = FALSE ~ 1 + fday, reverse = TRUE),
       data = fnc)
head(coef(fit2.ppom, mat = TRUE),  8)


# Miscellaneous output
head(fitted(fit2.ppom), 3)
head(predict(fit2.ppom), 3)
dim(model.matrix(fit2.ppom, type = "lm"))
dim(model.matrix(fit2.ppom, type = "vlm"))
constraints(fit2.ppom)[c(1,2,5,6)]


# ----------------------------------------------------------------------
# Marital Status Data
# Unfortunately the data set nzmarital is unavailable for distribution.


# ----------------------------------------------------------------------
# Back pain example.
# Requires the gnm package to be installed.

# install.packages("gnm")
data("backPain", package = "gnm")


set.seed(123)


# Scale the variables? Yes; the Anderson (1984) paper did (see his Table 6).
head(backPain, 4)
summary(backPain)
backPain <- transform(backPain,
                      sx1 = -scale(x1),
                      sx2 = -scale(x2),
                      sx3 = -scale(x3))


# Fit a rank-1 RR-MLM
set.seed(123)
bp.rrmlm1 <- rrvglm(pain ~ sx1 + sx2 + sx3, multinomial,
                    data = backPain)
Coef(bp.rrmlm1)  # The fitted A, C and B1 matrices
logLik(bp.rrmlm1)
}


# Fit a rank-2 RR-MLM (with a different normalization)
set.seed(123)
bp.rrmlm2 <- rrvglm(pain ~ sx1 + sx2 + sx3, multinomial, 
                    data = backPain, Rank = 2,
                    Corner = FALSE, Uncor = TRUE)

var(lv(bp.rrmlm2)) # Should be the same as diag(2)


# Biplot the data
 par(mfrow = c(1,1))
 par(mar = c(4.5,4.0,0.2,2.2)+0.1)
biplot(bp.rrmlm2, Acol = "blue", Ccol = "darkgreen", scores = TRUE,
       xlim = c(-4.5,2.2), ylim = c(-2.2, 2.2),
       chull = TRUE, clty = 2, ccol = "blue")


# ----------------------------------------------------------------------
# Miscellanea

iam(NA, NA, M = 4, both = TRUE, diag = TRUE)

# ----------------------------------------------------------------------