Bootstrap for a TVGEV Model

Source: R/TVGEV.R

Bootstrap for a TVGEV model.

# S3 method for class 'TVGEV'
bs(
  object,
  R = 100,
  type = c("param", "NP"),
  estim = "optim",
  parallel = FALSE,
  ...
)

Arguments

object

A TVGEV object.

R

Target number of bootstrap samples.

type

Character. The values "param" and "NP" can be used to chose between parametric and nonparametric bootstrap, see Details.

estim

Argument passed to MLE.

parallel

Logical. If TRUE the loop over the bootstrap sample is performed using foreach with %do par%.

...

Further arguments passed to MLE.

Value

A matrix containing the bootstraped parameter vectors, with one row for each bootstrap sample.

Details

The parametric bootstrap (type = "param") is as follows: for each bootstrap sample, a response vector \(\mathbf{y}^\star\) is drawn from the estimated GEV distribution as described in object; a MLE is performed to find the bootstraped parameter vector \(\boldsymbol{\psi}^\star\). For the nonparametric bootstrap (type = "NP"), the generalised residuals of the object as returned by residuals.TVGEV are resampled and are back-transformed to simulated observations; a MLE is then performed to find the bootstraped parameter vector \(\boldsymbol{\psi}^\star\) as in the parametric case.

Note

The simulated response does not cover the observations where the original response was NA.

Caution

For some bootstrap samples, the MLE may fail, in which case the coefficients will be ignored. So the true number of sampled coefficient vecors will generally be smaller than the target number as given in R.

Examples

example(TVGEV)
#> 
#> TVGEV> ## transform a numeric year into a date
#> TVGEV> df <- within(TXMax_Dijon, Date <- as.Date(sprintf("%4d-01-01", Year)))
#> 
#> TVGEV> df0 <- subset(df, !is.na(TXMax))
#> 
#> TVGEV> ## fit a TVGEV model. Only the location parameter is TV.
#> TVGEV> t1 <- system.time(
#> TVGEV+     res1 <- TVGEV(data = df, response = "TXMax", date = "Date",
#> TVGEV+                   design = breaksX(date = Date, breaks = "1970-01-01", degree = 1),
#> TVGEV+                   loc = ~ t1 + t1_1970))
#> 
#> TVGEV> ## The same using "nloptr" optimisation.
#> TVGEV> t2 <- system.time(
#> TVGEV+     res2 <- TVGEV(data = df, response = "TXMax", date = "Date",
#> TVGEV+                   design = breaksX(date = Date, breaks = "1970-01-01", degree = 1),
#> TVGEV+                   loc = ~ t1 + t1_1970,
#> TVGEV+                   estim = "nloptr",
#> TVGEV+                   parTrack = TRUE))
#> 
#> TVGEV> ## use extRemes::fevd the required variables need to be added to the data frame
#> TVGEV> ## passed as 'data' argument
#> TVGEV> t0 <- system.time({
#> TVGEV+    df0.evd <- cbind(df0, breaksX(date = df0$Date, breaks = "1970-01-01",
#> TVGEV+                     degree = 1));
#> TVGEV+    res0 <- fevd(x = df0.evd$TXMax, data = df0.evd, loc = ~ t1 + t1_1970)
#> TVGEV+  })
#> 
#> TVGEV> ## compare estimate and negative log-liks
#> TVGEV> cbind("fevd" = res0$results$par,
#> TVGEV+       "TVGEV_optim" = res1$estimate,
#> TVGEV+       "TVGEV_nloptr" = res2$estimate)
#>              fevd TVGEV_optim TVGEV_nloptr
#> mu0   32.06678895 32.06638460  32.06679282
#> mu1   -0.02391857 -0.02392656  -0.02391856
#> mu2    0.07727041  0.07728411   0.07727061
#> scale  1.75585289  1.75541862   1.75585353
#> shape -0.18130928 -0.18112018  -0.18131041
#> 
#> TVGEV> cbind("fevd" = res0$results$value,
#> TVGEV+       "VGEV_optim" = res1$negLogLik,
#> TVGEV+       "TVGEV_nloptr" = res2$negLogLik)
#>          fevd VGEV_optim TVGEV_nloptr
#> [1,] 177.2014   177.2014     177.2014
#> 
#> TVGEV> ## ====================================================================
#> TVGEV> ## use a loop on plausible break years. The fitted models
#> TVGEV> ## are stored within a list
#> TVGEV> ## ====================================================================
#> TVGEV> 
#> TVGEV> ## Not run: 
#> TVGEV> ##D 
#> TVGEV> ##D     yearBreaks <- c(1940, 1950, 1955, 1960:2000, 2005, 2010)
#> TVGEV> ##D     res <- list()
#> TVGEV> ##D 
#> TVGEV> ##D     for (ib in seq_along(yearBreaks)) {
#> TVGEV> ##D         d <- sprintf("%4d-01-01", yearBreaks[[ib]])
#> TVGEV> ##D         floc <- as.formula(sprintf("~ t1 + t1_%4d", yearBreaks[[ib]]))
#> TVGEV> ##D         res[[d]] <- TVGEV(data = df, response = "TXMax", date = "Date",
#> TVGEV> ##D         design = breaksX(date = Date, breaks = d, degree = 1),
#> TVGEV> ##D         loc = floc)
#> TVGEV> ##D     }
#> TVGEV> ##D 
#> TVGEV> ##D     ## [continuing...] ]find the model with maximum likelihood, and plot
#> TVGEV> ##D     ## something like a profile likelihood for the break date considered
#> TVGEV> ##D     ## as a new parameter. However, the model is not differentiable w.r.t.
#> TVGEV> ##D     ## the break! 
#> TVGEV> ##D 
#> TVGEV> ##D     ll <- sapply(res, logLik)
#> TVGEV> ##D     plot(yearBreaks, ll, type = "o", pch = 21, col = "orangered",
#> TVGEV> ##D          lwd = 2, bg = "gold", xlab = "break", ylab = "log-lik")
#> TVGEV> ##D     grid()
#> TVGEV> ##D     iMax <- which.max(ll)
#> TVGEV> ##D     abline(v = yearBreaks[iMax])
#> TVGEV> ##D     abline(h = ll[iMax] - c(0, qchisq(0.95, df = 1) /2),
#> TVGEV> ##D            col = "SpringGreen3", lwd = 2)
#> TVGEV> ##D 
#> TVGEV> ## End(Not run)
#> TVGEV> 
#> TVGEV> 
#> TVGEV> 
bsb <- bs(res2, R = 50, estim = "nloptr")
#> Warning: convergence not reached in optimisation
#> Warning: convergence not reached in optimisation

if (FALSE) { # \dontrun{
   library(parallel)
   library(doParallel)
   nc <- detectCores()
   cl <- makeCluster(nc)
   registerDoParallel(cl)

   ## findings: with 'nloptr', less than 1% of the optimisations
   ## diverges, while more than 10% diverged with 'optim'.
   te <- system.time(bspb <- bs(res2, R = 5000, estim = "nloptr",
                                parallel = TRUE))
   stopCluster(cl) 
} # }

## add the 'bs' results to the model object as 'boot' element,
## in order to facilitate the prediction. The number of
## bootstrap replicates in the prediction will correspond
## to the number of replications stored in the model object,
## an no longer to the default

res2$boot <- bsb
predict(res2, confintMethod = "boot")
#> Since 'object' is really time-varying, the Return Levels
#> depend on the date. A default choice of dates is made here.
#> Use the 'newdate' formal to change this.
#> 'object' embeds a bootstrap distribution. Formal arguments
#>  intented to be passed to `bs` (if any) will be ignored
#>          Date Period Level    Quant        L        U
#> 1  1920-01-01      5   95% 35.56868 34.40283 36.70862
#> 4  1920-01-01     10   95% 36.50728 35.10215 37.62304
#> 7  1920-01-01     20   95% 37.29520 35.64051 38.38083
#> 10 1920-01-01     50   95% 38.17369 36.19521 39.26103
#> 13 1920-01-01    100   95% 38.74128 36.52632 40.04032
#> 16 1920-01-01    150   95% 39.04056 36.69177 40.50191
#> 19 1920-01-01    200   95% 39.23966 36.79820 40.86563
#> 22 1920-01-01    500   95% 39.80798 37.08047 42.02426
#> 25 1920-01-01   1000   95% 40.17894 37.24078 42.80501
#> 2  1970-01-01      5   95% 34.37272 33.61482 35.08380
#> 5  1970-01-01     10   95% 35.31131 34.43236 36.12135
#> 8  1970-01-01     20   95% 36.09924 34.96320 36.88258
#> 11 1970-01-01     50   95% 36.97773 35.57076 38.05899
#> 14 1970-01-01    100   95% 37.54532 35.87252 38.84859
#> 17 1970-01-01    150   95% 37.84460 36.01994 39.24527
#> 20 1970-01-01    200   95% 38.04370 36.11345 39.52152
#> 23 1970-01-01    500   95% 38.61202 36.38916 40.46518
#> 26 1970-01-01   1000   95% 38.98298 36.55484 41.27275
#> 3  2020-01-01      5   95% 37.04025 35.99741 37.99838
#> 6  2020-01-01     10   95% 37.97884 36.80525 38.98024
#> 9  2020-01-01     20   95% 38.76677 37.45658 39.85426
#> 12 2020-01-01     50   95% 39.64526 38.18984 40.97480
#> 15 2020-01-01    100   95% 40.21285 38.42933 41.79163
#> 18 2020-01-01    150   95% 40.51213 38.53748 42.23067
#> 21 2020-01-01    200   95% 40.71123 38.61801 42.53170
#> 24 2020-01-01    500   95% 41.27955 38.90997 43.43640
#> 27 2020-01-01   1000   95% 41.65051 39.08470 44.06985