Nothing
plotrl.evmOpt <- # intended as a diagnostic for an evm fitted with no covariates. Called by plot.evm
function(object, alpha = .050,
xlab, ylab, main,
pch= 1, ptcol =2 , cex=.75, linecol = 4 ,
cicol = 0, polycol = 15, smooth = FALSE, RetPeriodRange=NULL ){
wh <- sapply(object$data$D, ncol)
if (any(wh > 1)){
stop("use plot method for object returned by predict.evmOpt to see rl plots if covariates in model")
}
if (missing(xlab) || is.null(xlab)) { xlab <- "Return period" }
if (missing(ylab) || is.null(ylab)) { ylab <- "Return level" }
if (missing(main) || is.null(main)) { main <- "Return Level Plot" }
wh <- getPlotRLdata(object, alpha, RetPeriodRange)
plotRLevm(wh$m[wh$plotX], wh$xm[wh$plotX,],
polycol, cicol, linecol, ptcol,
wh$n, wh$xdat, pch, smooth, xlab, ylab, main,
xrange=wh$xrange, yrange=wh$yrange)
invisible(list(m=wh$m, xm=wh$xm))
}
getPlotRLdata <- function(object, alpha, RetPeriodRange){
xdat <- object$data$y
n <- length(xdat) / object$rate # Number of obs prior to thresholding
if(!is.null(RetPeriodRange)){
ran <- log10(RetPeriodRange)
jj <- seq(ran[1],ran[2],by=0.1)
} else {
jj <- seq(log10(1/(1-1/n)), 2*log10(n),by=0.1)
}
tol <- 0.00001 # to avoid div by zero in gev rl calc
m <- unique(c(max(1/object$rate,1+tol), 10^jj)) / object$rate
xm <- matrix(unlist(rl(object, M=m, ci.fit=TRUE, alpha=alpha)), ncol=3, byrow=TRUE)
U <- object$threshold - abs(object$threshold/100)
plotX <- xm[,1] > U
xrange <- range(m)
yrange <- range(c(xdat, range(xm[plotX,])))
list(n=n, xdat=xdat, m=m, xm=xm, plotX=plotX, xrange=xrange, yrange=yrange)
}
#' @rdname rl
#' @param plot. Parameter for plot method, whether to produce plots.
#' @export
plot.rl.evmOpt <- function(x, # method for rl.(evmBoot or evmSim or evmOpt) object, which may have covariates. Plots return level for each unique row in design matrix
xlab, ylab, main,
pch= 1, ptcol =2 , cex=.75, linecol = 4 ,
cicol = 0, polycol = 15, smooth = FALSE, sameAxes=TRUE, type="median", ylim = NULL, plot.=TRUE, ...){
Class <- class(x)
x <- x$obj
if (missing(xlab) || is.null(xlab)) { xlab <- "Return period" }
if (missing(ylab) || is.null(ylab)) { ylab <- "Return level" }
if (missing(main) || is.null(main)) {
main <- "Return Level Plot"
SetMain <- TRUE
} else {
if(length(main) == 1){
SetMain <- TRUE
} else {
SetMain <- FALSE
}
}
nm <- length(names(x))
if(nm < 2) {
stop("Need to have more than one value of M at which to plot return level curve")
}
nd <- dim(x[[1]])[2]
ncov <- length(unlist(x)) / (nm * nd)
ValNames <- colnames(x[[1]])
if(length(ValNames) == 1 | length(grep("%",ValNames)) < 2){
stop("Please use ci.fit=TRUE in call to predict, to calculate confidence intervals")
}
if(!SetMain & length(main) != ncov){
stop("main must be length 1 or number of unique covariates for prediction")
}
if(Class == "rl.evmOpt"){
if(any(colnames(x[[1]]) == "se.fit")){
which <- colnames(x[[1]]) != "se.fit"
nd <- nd-1
ValNames <- ValNames[which]
Unlist <- unlist(x)[rep(which,each=ncov)]
} else {
Unlist <- unlist(x)
}
} else if(Class == "rl.evmSim" |Class == "rl.evmBoot"){
if(casefold(type) == "median"){
which <- substring(colnames(x[[1]]), nchar(colnames(x[[1]])) -3) != "Mean"
} else if(casefold(type) == "mean") {
which <- substring(colnames(x[[1]]), nchar(colnames(x[[1]])) -2) != "50%"
} else {
stop("type must be \"mean\" or \"median\" ")
}
nd <- nd-1
ValNames <- ValNames[which]
Unlist <- unlist(x)[rep(which,each=ncov)]
}
Array <- array(Unlist,c(ncov,nd,nm),dimnames=list(NULL,ValNames,names(x)))
m <- as.numeric(substring(dimnames(Array)[[3]],first=3))
if(ncov>1){
covnames <- dimnames(Array)[[2]][-(1:3)]
} else {
covnames <- ""
}
if(sameAxes){
if(!is.null(ylim)){
yrange <- ylim
} else {
yrange <- range(Array[,1:3,])#ifelse(is.null(ylim),,ylim)
}
}
for(i in 1:ncov){
xm <- t(Array[i,,])
cov <- xm[1,-(1:3)]
if(!sameAxes){
yrange <- range(xm[,1:3])
}
if(SetMain){
if(length(covnames) == 1 && covnames != ""){
Main <- paste(main,"\n", paste(covnames,"=",signif(cov,4),collapse=", "))
} else {
Main <- main
}
} else {
Main <- main[i]
}
if(plot.){
plotRLevm(m,xm,polycol = polycol,cicol=cicol,linecol=linecol,ptcol=ptcol,pch=pch,
smooth=smooth,xlab=xlab,ylab=ylab,main=Main,xrange=range(m),yrange=yrange)
}
}
invisible(list(m=m,xm=Array))
}
#' @rdname rl
#' @export
plot.rl.evmSim <- plot.rl.evmOpt
#' @rdname rl
#' @export
plot.rl.evmBoot <- plot.rl.evmOpt
plotRLevm <- function(M,xm,polycol,cicol,linecol,ptcol,n,xdat,pch,smooth,xlab,ylab,main,xrange,yrange){
# worker function - called by plotrl.evmOpt, plot.rl.evmOpt, plot.rl.evmSim, plot.rl.evmBoot
o <- order(M) # in case the return period are not in ascending order.
M <- M[o]
xm <- xm[o,]
plot(M, xm[,1], log = "x", type = "n",
xlim=xrange, ylim=yrange, xlab = xlab, ylab = ylab, main = main)
if (smooth) {
splo <- spline(log(M), xm[,2], 200)
sphi <- spline(log(M), xm[,3], 200)
if ( polycol != 0 ) {
polygon( exp(c( splo$x, rev( sphi$x ) )),
c( splo$y, rev( sphi$y ) ),
col = polycol , border=FALSE)
} # Close if (polycol
lines( exp(splo$x), splo$y, col = cicol )
lines( exp(sphi$x), sphi$y, col = cicol )
} else{
if (polycol != 0){
polygon(c(M, rev( M)),
c(xm[,2],rev(xm[,3])),
col=polycol, border = FALSE) # Close polygon
} else {
lines(M, xm[,2], col = cicol)
lines(M, xm[,3], col = cicol)
}
}
lines(M, xm[,1], col = linecol[ 1 ] )
# Add observed data to the plot
if(!missing(xdat) & !missing(n)){
ly <- length(xdat)
points(1 / (1 - ((n - ly + 1):n) / (n + 1)), sort(xdat), pch=pch, col=ptcol)
box()
}
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.