# Analysis of the 2008 WFFC data
# T. W. Yee (2013)
# Fit a bivariate (logistic) odds-ratio model (BOM).

# ======================================================================
# Read the raw data in

library("VGAM")      # Requires Version 0.9-2 or later
library("VGAMdata")  # Requires Version 0.9-2 or later



data(wffc)
summary(wffc)
dim(wffc)



# ======================================================================
# Fit a BOM. 
# Y1 = Waihou (1) vs. Waimakariri (0)
# Y2 = mornaft = morning (0) vs. afternoon (1)
# Log odds ratio is significantly different from 0.


# Process Waihou River
waihou <- subset(wffc, water == "Waihou")
temp.waihou <- na.omit(waihou)
temp.waihou <- subset(temp.waihou,
                      session %in% c(1,2,5,6)) # Days 1 or 3 only
temp.waihou <- transform(temp.waihou,
                         day1    = ifelse(session %in% 1:2, 1, 0),
                         mornaft = 1 - session %% 2)
summary(temp.waihou)


# Process Waimakariri River
waimak <- subset(wffc, water == "Waimakariri")
temp.waimak <- na.omit(waimak)

temp.waimak <- subset(temp.waimak,
                      session %in% c(1,2,5,6)) # Days 1 or 3 only
temp.waimak <- transform(temp.waimak,
                         day1    = ifelse(session %in% 1:2, 1, 0),
                         mornaft = 1 - session %% 2)
summary(temp.waimak)



# Merge the two rivers into a data frame
waiwai <- rbind(temp.waihou, temp.waimak)
waiwai <- transform(waiwai,
                    loglength = log(length),
                    lengthcm  = length/10)


# Fit the model
fit.blom <- vglm(cbind(water == "Waihou", mornaft) ~ ns(lengthcm, df = 3),
                 binom2.or(zero = 3), data = waiwai, trace = TRUE, crit = "c")


# Plot the model
par(mfrow = c(1,2), las = 1, oma = c(1, 0, 0, 0), cex = 1.2)
plotvgam(fit.blom, se = T, scol = "blue", scale = 4, xlab = "Length", which.cf = 1)
text(x = 20, y = 2.4, "(a)", font = 1)

plotvgam(fit.blom, se = T, scol = "blue", scale = 4, xlab = "Length", which.cf = 2)
text(x = 20, y = 3.17, "(b)", font = 1)



# Look at some results 
colSums(fit.blom@y)
(sfit.blom <- summary(fit.blom))   # oratio may well be 1
coef(fit.blom, matrix = TRUE)


# Estimated odds ratio
exp( coef(fit.blom, matrix = TRUE)["(Intercept)","log(oratio)"] )

2*pnorm(-abs(sfit.blom@coef3["(Intercept):3", "t value"])) # 2-sided pvalue

# ======================================================================