# Read the rainfall data into "ysrain"
source("http://www.stat.auckland.ac.nz/~yee/updates/ysrain.R")

# Alternatively, download the file and then run the following
#if( .Platform$OS.type == "windows") {
#    source(file.choose()) # Windows; choose ysrain.R
#} else {
#    source("./ysrain.R") # Unix or Linux
#}



library("VGAM")  # Version 0.9-4 or later is required





threshold <- ysrain$Thresh[ysrain$CluMax]
clmax <- ysrain$CluMax
rainfit <-
  vglm(CenSite ~ bs(SiteA) + bs(SiteB) + bs(SiteC) + bs(SiteD) + Year,
       gpd(threshold, lshape = "identitylink", zero = 2),
       data = ysrain, subset = clmax, maxit = 50)

par(mfrow = c(3, 2), las = 1, mar = c(5, 4, 1, 1)+0.1)
plotvgam(rainfit, se = TRUE, lcol = "blue", scol = "darkgreen",
         ylim = c(-0.5, 0.8))

coef(rainfit, matrix = TRUE)

exp(8 * coef(rainfit)["Year"])  # Approx 0.9

summary(rainfit)