# fremantle.R
# Example code from Yee and Stephenson (2007), Extremes.
# Last modified: 2018-11-05

# Assign "fremantle" the data frame by the same name in the ismev package.
# Note that old versions of ismev had the variable name "Sea Level" rather
# than "SeaLevel" and this causes problems when used within a formula.
# library("ismev"); data("fremantle"); fremantle = fremantle; detach()

data("fremantle", package = "ismev")




# The latest version is at http://www.stat.auckland.ac.nz/~yee/VGAM
library("VGAM")


# Model 1
fit1 <- vgam(SeaLevel ~ s(Year, df=3) + s(SOI, df=3), 
            family=gev(zero="shape"),             # New
#           family=gev(lshape="logoff", zero=3),  # Old
            data=fremantle, maxit=50)
par(mar=c(5,4.3,2.5,1.3)+0.1, mfrow=c(2,2), las=1)
plot(fit1, se=TRUE, lcol="blue", scol="limegreen")

summary(fit1)  # 3 of the p-values for linearity are > 0.20
coef(fit1, matrix=TRUE)
logoff(-1.493495, offset=0.5, inverse = TRUE)  # Estimate of xi
#logoff(-1.493495, list(offset=0.5), inv=TRUE)  # Estimate of xi



# Model 2
cm <- list("(Intercept)"=diag(3), "Year"=rbind(1,0,0),
          "bs(Year, degree = 1, knot = 1945)"=rbind(0,1,0),
          "SOI"=diag(3)[,1:2])
fit2 <- vglm(SeaLevel ~ Year + bs(Year, degree=1, knot=1945) + SOI,
            family=gev(zero="shape"),             # New
#           family=gev(lshape="logoff", zero=3),  # Old
            data=fremantle, constraints=cm)
par(mfrow=c(2,2), mar=c(5,4,2.5,0.1)+0.1, las=1)
plot(as(fit2, "vgam"), se=TRUE, lcol="blue", scol="limegreen")
coef(fit2, matrix=TRUE)
1 - pchisq(2*(logLik(fit1)-logLik(fit2)), 
           df=df.residual(fit2)-df.residual(fit1)) # Deviance test