Description Details Author(s) References Examples
Conditional multivariate extreme values modelling following the approach of Heffernan and Tawn, 2004.
Package: | mex |
Type: | Package |
Version: | 1.2 |
Date: | 2011-10-06 |
License: | GPL (>=2) | BSD |
gpd
: Fit generalized Pareto distributions to data, possibly with
covariates. Use maximum likelihood estimation, maximum
penalized likelihood estimation or simulate from the posterior
distribution.
mex
: Fit multiple, independent generalized Pareto models
to the the upper tails of the columns of a data set, and estimate the conditional
dependence structure between the columns using the method of Heffernan and
Tawn.
bootmex
: Bootstrap estimation for parameters in generalized Pareto models
and in the dependence structure.
declust
: Estimation of extremal index and subsequent declustering of dependent sequences using the intervals estimator of Ferro and Segers.
validate.texmex
: Use the RUnit
package to run lots of tests to
establish a high degree of confidence that the functions do what
they are supposed to.
Harry Southworth, Janet E. Heffernan
Maintainer: Harry Southworth <harry.southworth@astrazeneca.com>
J. E. Heffernan and J. A. Tawn, A conditional approach for multivariate extreme values, Journal of the Royal Statistical society B, 66, 497 – 546, 2004.
C.A.T Ferro and J. Segers, Inference for Clusters of Extreme Values, Journal of the Royal Statistical society B, 65, 545 – 556, 2003.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 | # Analyse the winter data used by Heffernan and Tawn
mymex <- mex(winter, mqu = .7, penalty="none", dqu=.7, which = "NO")
plot(mymex)
myboot <- bootmex(mymex)
plot(myboot)
mypred <- predict(myboot, pqu=.95)
summary(mypred , probs = c( .025, .5, .975 ))
# Analyse the liver data included in the package
library(MASS) # For the rlm function
liver <- liver[liver$ALP.M > 1,] # Get rid of outlier
liver$ndose <- as.numeric(liver$dose)
alt <- resid(rlm(log(ALT.M) ~ log(ALT.B) + ndose, data=liver, method="MM"))
ast <- resid(rlm(log(AST.M) ~ log(AST.B) + ndose, data=liver, method="MM"))
alp <- resid(rlm(log(ALP.M) ~ log(ALP.B) + ndose, data=liver, method="MM"))
tbl <- resid(rlm(log(TBL.M) ~ log(TBL.B) + ndose, data=liver, method="MM"))
r <- data.frame(alt=alt, ast=ast, alp=alp, tbl=tbl)
Amex <- mex(r[liver$dose == "A",], mqu=.7)
Bmex <- mex(r[liver$dose == "B",], mqu=.7)
Cmex <- mex(r[liver$dose == "C",], mqu=.7)
Dmex <- mex(r[liver$dose == "D",], mqu=.7)
par(mfcol=c(3,3))
plot(Amex)
plot(Dmex, col="blue")
## Take a closer look at the marginal behaviour of ALT
r$ndose <- liver$ndose
altmod1 <- gpd(alt, qu=.7, phi = ~ ndose, xi = ~ ndose, data=r)
altmod2 <- gpd(alt, qu=.7, phi = ~ ndose, data=r)
altmod3 <- gpd(alt, qu=.7, xi = ~ ndose, data=r)
altmod4 <- gpd(alt, qu=.7, data=r)
# Prefer model 3, with term for xi on basis of AIC
# Lines commented out to keep CRAN robots happy
#balt3 <- gpd(alt, qu=.7, xi = ~ ndose, data=r, method="simulate")
#par(mfrow=c(3,3))
#plot(balt3)
# use longer burn-in and also thin the output
#balt3 <- thinAndBurn(balt3,burn=1000,thin=5)
#plot(balt3)
# Get some simulated values for dose D
#DParam <- predict(balt3,type="lp",newdata=data.frame(ndose=4),all=TRUE)[[1]]
#simD <- rgpd(nrow(DParam), sigma=exp(DParam[,"phi"]), xi=DParam[,"xi"], u=quantile(alt, .7))
# These are simulated residuals. Get some baselines and transform all
# to raw scale
#b <- sample(log(liver$ALT.M), size=nrow(balt3$param), replace=TRUE)
#res <- exp(b + simD)
# estimate quantiles on raw scale
#quantile(res, prob=c(.5, .75, .9, .95, .99))
# estimate proportion exceeding 3*upper limit of normal
#mean(res > 36 * 3) # 36 is the upper limit of normal for ALT
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.