Venice example and data
This example was studied by Coles and Fawcett who both fitted time varying GEV models with linear and quadratic trend by using the ismev package. We reproduce here their results by using NSGEV and also extRemes.
Although the data contains the \(r\)-largest observations, we use only the annual maxima here; These are found in column 2 of the matrix object.
data(venice, package = "ismev")
df <- data.frame(Year = venice[ , 1],
date = as.Date(sprintf( "%4d-01-01", venice[ , 1])),
SeaLevel = venice[ , 2])
df <- within(df, date <- as.Date(sprintf("%4d-01-01", Year)))
plot(SeaLevel ~ date, data = df, type = "o", pch = 16, xlab = "")
Fitting stationary GEV models
We begin by fitting a stationary GEV model.
estim <- "optim"
fitTVGEV0 <- TVGEV(data = df,
response = "SeaLevel", date = "date",
estim = estim,
loc = ~ 1)
fitismev0 <- gev.fit(xdat = df$SeaLev)
fitextRemes0 <- extRemes::fevd(x = df$SeaLev, type = "GEV")The estimated coefficients and maximised log-likelihoods are as follows.
## mu_0 sigma_0 xi_0
## NSGEV 111.0982 17.17421 -0.07672147
## ismev 111.0993 17.17549 -0.07673265
## extRemes 111.0973 17.17910 -0.07678211## NSGEV ismev extRemes
## -222.7145 -222.7145 -222.7145Linear trend
We now introduce a linear trend. We first create a dataframe to be used with the functions from extRemes and ismev, with a linear and variable a quadratic trend variable as well.
dfismev <- within(df, {t1num <- Year - 1930; t2num <- (Year - 1930)^2 })
fitTVGEV1 <- TVGEV(data = df,
response = "SeaLevel", date = "date",
design = polynomX(date = date,
degree = 2, origin = "1930-01-01"),
estim = estim,
trace = 2,
loc = ~ 1 + t1)
fitismev1 <- gev.fit(xdat = df$SeaLev,
ydat = as.matrix(dfismev[ , "t1num"]),
mul = 1)
fitextRemes1 <- extRemes::fevd(x = dfismev$SeaLev,
data = dfismev,
type = "GEV",
location.fun = ~t1num)## mu_0 mu_t1 sigma_0 xi_0
## NSGEV 96.97787 0.5645199 14.58559 -0.02750053
## ismev 96.98579 0.5641427 14.58435 -0.02731421
## extRemes 96.96317 0.5648220 14.58120 -0.02731842## NSGEV ismev extRemes
## -216.0625 -216.0626 -216.0626We can perform a likelihood-ratio test as in the references.
anova(fitTVGEV0, fitTVGEV1)## Analysis of Deviance Table
##
## df deviance W Pr(>W)
## fitTVGEV0 3 445.43
## fitTVGEV1 4 432.13 13.304 0.0002648 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Quadratic trend
We now add a quadratic term in the trend.
fitTVGEV2 <- TVGEV(data = df,
response = "SeaLevel", date = "date",
design = polynomX(date = date,
degree = 2, origin = "1930-01-01"),
estim = estim,
trace = 2,
loc = ~ 1 + t1 + t2)
fitismev2 <- gev.fit(xdat = dfismev$SeaLev,
ydat = as.matrix(dfismev[ , c("t1num", "t2num")]),
mul = c(1, 2))
fitextRemes2 <- extRemes::fevd(x = dfismev$SeaLev,
data = dfismev,
type = "GEV",
location.fun = ~t1num + t2num) The estimated coefficients and maximised log-likelihoods are as follows.
## mu_0 mu_t1 mu_t2 sigma_0 xi_0
## NSGEV 98.23632 0.4398647 0.002236958 14.62472 -0.02896788
## ismev 96.38501 0.6325000 -0.001340599 14.56492 -0.02547159
## extRemes 96.37645 0.6330981 -0.001347965 14.56132 -0.02545380## NSGEV ismev extRemes
## -216.1088 -216.0555 -216.0555So, compared to the two other packages, NSGEV::TVGEV
returns different estimated coefficients and a smaller maximised
log-likelihood, indicating a problem in the maximisation. Since the
scaling of the trend variable is often a concern, we try a change in the
time origin.
fitTVGEV2b <- TVGEV(data = df,
response = "SeaLevel", date = "date",
design = polynomX(date = date,
degree = 2, origin = "1950-01-01"),
estim = estim,
trace = 2,
loc = ~ 1 + t1 + t2)
rbind("NSGEV" = coef(fitTVGEV2b),
"ismev" = fitismev2$mle,
"extRemes" = fitextRemes2$results$par)## mu_0 mu_t1 mu_t2 sigma_0 xi_0
## NSGEV 108.50613 0.5793764 -0.001374096 14.55759 -0.02535797
## ismev 96.38501 0.6325000 -0.001340599 14.56492 -0.02547159
## extRemes 96.37645 0.6330981 -0.001347965 14.56132 -0.02545380
c("NSGEV" = logLik(fitTVGEV2b),
"ismev" = -fitismev2$nllh,
"extRemes" = -fitextRemes2$results$value)## NSGEV ismev extRemes
## -216.0554 -216.0555 -216.0555The log-likelihood maximised with NSGEV::TVGEV is now
identical to that of the other two packages so a simple change of origin
solved the problem. The resulting model in the TVGEV object
is simply a re-parameterisation of the former, and the comparison with
the other packages needs some care because the trends for the location
parameter are no longer identical. The new parameterisation facilitates
the maximisation. To complete the comparison we can use 1950 as an
origin for all models.
dfismevb <- within(df, {t1num <- Year - 1950; t2num <- (Year - 1950)^2 })
fitismev2b <- gev.fit(xdat = dfismevb$SeaLev,
ydat = as.matrix(dfismevb[ , c("t1num", "t2num")]),
mul = c(1, 2))
fitextRemes2b <- extRemes::fevd(x = dfismevb$SeaLev,
data = dfismevb,
type = "GEV",
location.fun = ~t1num + t2num)## mu_0 mu_t1 mu_t2 sigma_0 xi_0
## NSGEV 108.5061 0.5793764 -0.001374096 14.55759 -0.02535797
## ismev 108.5073 0.5798072 -0.001392904 14.55171 -0.02536854
## extRemes 108.5062 0.5791653 -0.001358444 14.56185 -0.02550671
c("NSGEV" = logLik(fitTVGEV2b),
"ismev" = -fitismev2b$nllh,
"extRemes" = -fitextRemes2b$results$value)## NSGEV ismev extRemes
## -216.0554 -216.0555 -216.0555We now can proceed to the likelihood ratio test.
anova(fitTVGEV1, fitTVGEV2b)## Analysis of Deviance Table
##
## df deviance W Pr(>W)
## fitTVGEV1 4 432.13
## fitTVGEV2b 5 432.11 0.014213 0.9051The test tells that adding a quadratic trend term does not improve the model.
Findings
With this example it appears that in NSGEV::TVGEV the
maximisation of log-likelihood can fail to find the best value, still
giving an indication of convergence. Here the problem is due to the
order of magnitude of one covariate, namely the quadratic term in the
polynomial basis. It was fixed by simply changing the origin of
time.
With the current implementation of NSGEV when formulas involving several basis covariates are used, we recommend to fit models with increasing complexity, and check that the log-likelihood increases when adding a new covariate.
Several solutions can be to improve the code of
NSGEV::TVGEV in future versions.
Allow the user to give initial values for the parameters.
Try a “multistart” optimisation with randomly selected initial values.
Scale the design matrix \(\mathbf{X}\) by subtracting its mean to each column, and by dividing each column by its standard deviation. Although this is a quite simple application of the
scalemethod for the"matrix"class, the estimated coefficients should then be translated back in the original scale, which can be quite tedious to implement.