options(prompt = "> ", continue="  ", useFancyQuotes = FALSE)
options(width = 80)
options(digits = 5)




## ---------------------------------------------------------------------------
library("VGAM")
library("VGAMdata")
 ps.options(pointsize = 12, family = "Times")
pdf.options(pointsize = 12, family = "Times")







## ------------------------------------------------------------------------
# Note: this is computationally expensive
set.seed(1)
hspider[, 1:6] <- scale(hspider[, 1:6])  # Standardize the environmental variables
p1cao.hs <- cao(cbind(Alopacce, Alopcune, Alopfabr, Arctlute, Arctperi, Auloalbi,
                      Pardlugu, Pardmont, Pardnigr, Pardpull, Trocterr, Zoraspin) ~
                WaterCon + BareSand + FallTwig + CoveMoss + CoveHerb + ReflLux,
                poissonff, data = hspider,
                trace = FALSE,  # The difference with the above
                df1.nl = 2, Bestof = 10, Crow1positive = FALSE)

## ------------------------------------------------------------------------
par(mar = c(3.4, 3, 1.4, 0.2) + 0.1, mfrow = c(4, 3))
par(mgp = c(3, 1, 0))  # Default: c(3, 1, 0)

plot(p1cao.hs, lcol = "blue", lwd = 2, ylim = c(-5, 5), xlab = "", ylab = "")

## ------------------------------------------------------------------------
round(concoef(p1cao.hs), digits = 2)









## ------------------------------------------------------------------------
# Note: this is computationally expensive
trapO <- transform(trapO, sc.doy     = scale(doy),
                          sc.LevelTW = scale(LevelTW),
                          sc.MinAT   = scale(MinAT),
                          sc.MaxAT   = scale(MaxAT),
                          f.year     = factor(Year))
set.seed(123)  # For reproducibility
p1cao.to <- cao(cbind(BFTW, BMTW, RFTW, RMTW) ~ 
                sc.doy + f.Year + sc.MinAT + sc.MaxAT + sc.LevelTW,
                family = poissonff, df1.nl = 2.5, trace = FALSE, data = trapO)


Coef(p1cao.to)


par(mar = c(3.4, 3, 1.4, 0.2) + 0.1, mfrow = c(1, 1))
index <- 1:ncol(depvar(p1cao.to))
persp(p1cao.to, col = index, label = TRUE, las = 1)







## ------------------------------------------------------------------------
# Note: this is computationally expensive
xs.nz.em <- subset(xs.nz, 45 < age  & age < 55 & sex == "M" & smokenow == 0 &
                  ethnicity == "European" & fh.cancer == 0 & fh.heartdisease == 0)
sort(colMeans(xs.nz.em[, c("asthma","cancer","diabetes","heartattack","stroke")],
              na.rm = TRUE))  # Disease prevalences

set.seed(123)
b1cao.xs <-
  cao(cbind(asthma, cancer, diabetes, heartattack, stroke) ~
      depressed + embarrassed + fedup + hurt + miserable +   # 11 psychological
      nofriend + moody + nervous + tense + worry + worrier,  # variables
      df1.nl = 1.25, Crow1positive = FALSE, trace = FALSE,
      binomialff(multiple.responses = TRUE), data = xs.nz.em)
nrow(depvar(b1cao.xs))  # n
sort(deviance(b1cao.xs, history = TRUE))


round(sort(concoef(b1cao.xs)[, 1]), digits = 2)


par(mar = c(5.1, 4, 1.7, 1.5) + 0.1)
par(mfrow = c(2, 3))
plot(b1cao.xs, ylim = c(-10, 0), lcol = "blue", las = 1)


lvdata <- data.frame(latvar1 = c(latvar(b1cao.xs)))
check.b1cao.xs <- 
  vgam(depvar(b1cao.xs) ~ s(latvar1, df = 1.25), # Must be same as the original
       binomialff(multiple.responses = TRUE), data = lvdata, trace = FALSE)
par(mfrow = c(2, 5))
plot(check.b1cao.xs, lcol = "blue", scol = "orange", se = TRUE, ylim = c(-4, 4),
     las = 1)
plot(check.b1cao.xs, lcol = "blue", deriv = 1)