Confidence intervals for a TVGEV object
An object with class "TVGEV".
Parameter name. NOT USED YET.
Confidence level.
"delta", "proflik",
"boot".
Integer level of verbosity.
Logical. If TRUE the confidence limits will be
rounded to a small number of digits. This number is chosen using the
smallest of the standard deviations for the estimated parameters.
Character giving the class of the output. By default a
three-dimensional array with dimensions: parameter,
lower/upper limit, and level. If level has length
1, using drop on the output will remove the third dimension
and produce a matrix. The "data.frame" gives the same results
in 'long' format.
Arguments to be passed to the profLik or to the
bs method. For instance, they can be used to chose the type
of bootstrap or the number of bootstrap replications if
method is "boot".
An array or a data frame with the lower and upper bounds
"L" and "U" of the confidence intervals.
For the bootstrap method(s), the time required depends on
whether object embeds a bootstrap distribution or not. When
no bootstrap distribution is found, it has to be computed using
the bs method; the formal arguments that are given in
... and intended to be used by bs will then be
ignored with a warning.
## use the example of TVGEV which defines a TVGEV object named
## 'res2'
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.06679233
#> mu1 -0.02391857 -0.02392656 -0.02391860
#> mu2 0.07727041 0.07728411 0.07727031
#> scale 1.75585289 1.75541862 1.75585346
#> shape -0.18130928 -0.18112018 -0.18130938
#>
#> 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>
## (default) delta method, several confidence levels
ci <- confint(res2, level = c(0.90, 0.70))
ci1 <- confint(res2, level = c(0.95, 0.90, 0.70), out = "data.frame")
ci2 <- confint(res2, level = 0.95, out = "data.frame")
## Profile-likelihood
ci3 <- confint(res2, level = 0.95, method = "proflik")
#>
#>
#> o Finding CI for "mu_0"
#>
#> o 95%, lower bound: 31.28
#> Constraint check -0.0000002, 179.1221
#> 95%, upper bound: 32.84
#> Constraint & value -0.0000001, 179.1221
#>
#>
#> o Finding CI for "mu_t1"
#>
#> o 95%, lower bound: -0.02
#> Constraint check -1.9207294, 177.2014
#> 95%, upper bound: -0.00
#> Constraint & value -0.0008238, 179.1213
#>
#>
#> o Finding CI for "mu_t1_1970"
#>
#> o 95%, lower bound: 0.05
#> Constraint check -0.0243793, 179.0977
#> 95%, upper bound: 0.10
#> Constraint & value -0.0001081, 179.1220
#>
#>
#> o Finding CI for "sigma_0"
#>
#> o 95%, lower bound: 1.50
#> Constraint check -0.0000000, 179.1221
#> 95%, upper bound: 2.09
#> Constraint & value -0.0000001, 179.1221
#>
#>
#> o Finding CI for "xi_0"
#>
#> o 95%, lower bound: -0.31
#> Constraint check -0.0000008, 179.1221
#> 95%, upper bound: -0.03
#> Constraint & value -0.0000000, 179.1221