R/plot.R
In philmu/climex: Estimating Extreme Events in Climatic Time Series (using Data of the DWD)
Defines functions
plot.climex.fit.gpd
plot.climex.fit.gev
likelihood.plot
ttplot
multiplot
Documented in
likelihood.plot
multiplot
plot.climex.fit.gev
plot.climex.fit.gpd
ttplot
### plot.R - Functions for plotting either the 'xts' class data or the
## fitting results
##' @title Plot multiple time series.
##' @description This function takes several \pkg{ggplot2} plots and
##' arranges them on one page.
##' @details Well, it is not really related to extreme value fitting
##' and I am aware of the fact that its quite bad style to export
##' auxiliary functions. But this one is just so handy. If the
##' \strong{layout} is something like
##' \code{matrix( c( 1, 2, 3, 3 ), nrow = 2, byrow = TRUE )},
##' then plot 1 will go in the upper left, 2 will go in the upper
##' right, and 3 will go all the way across the bottom. ggplot
##' objects can be passed in ..., or to \strong{plotlist} (as a list
## of \emph{ggplot} objects).
##'
##' @param plotlist Alternative way to pass \emph{ggplot} objects.
##' @param tt.title Same as \strong{main}.
##' @param main Title of the overall plot.
##' @param ... Additional arguments for ggplot.
##' @param cols Number of columns in layout. Default = 1.
##' @param layout A matrix specifying the layout. If present,
##' \strong{cols} is ignored.
##'
##' @family plot
##' @author Paul Teetor (original), Philpp Mueller (minor changes)
##' @examples
##' anomalies.potsdam <- anomalies( temp.potsdam )
##' multiplot( anomalies.potsdam, temp.potsdam, cols = 2, main =
##' "Difference between the pure data of the daily maximum
##' temperatures of the Potsdam station and their anomalies" )
##' @return print object.
##' @export
##' @import ggplot2
##' @import grid
multiplot <- function( plotlist = NULL, tt.title = main, main = NULL,
..., cols = 1, layout = NULL ) {
## Make a list from the ... arguments and plotlist
plots <- c( list( ... ), plotlist )
numPlots = length( plots )
## If layout is NULL, then use 'cols' to determine layout
if ( is.null( layout ) ) {
## Make the panel
## ncol: Number of columns of plots
## nrow: Number of rows needed, calculated from # of cols
layout <- matrix( seq( 1, cols* ceiling ( numPlots/ cols ) ),
ncol = cols, nrow = ceiling( numPlots/ cols ) )
}
if ( numPlots == 1 ) {
print( plots[[ 1 ]] )
} else {
## Set up the page
grid.newpage()
if ( !is.null( tt.title ) ){
pushViewport( viewport(
layout = grid.layout(
nrow( layout ) + 1, ncol( layout ),
heights = unit( rep( 1, nrow( layout ) + 1 ),
c( "lines", rep( "null",
nrow( layout) ) ) ) ) ) )
} else
pushViewport( viewport(
layout = grid.layout( nrow( layout ), ncol( layout ) ) ) )
## Make each plot, in the correct location
if ( !is.null( tt.title ) ){
grid.text( tt.title, vp = viewport(
layout.pos.row = 1,
layout.pos.col = c( 1, ncol( layout ) )
) )
for ( ii in 1 : numPlots ) {
## Get the i,j matrix positions of the regions that contain
## this subplot
matchidx <- as.data.frame( which( layout == ii,
arr.ind = TRUE ) )
print( plots[[ ii ]], vp = viewport(
layout.pos.row = matchidx$row + 1,
layout.pos.col = matchidx$col ) )
}
} else {
for ( ii in 1 : numPlots ) {
## Get the i,j matrix positions of the regions that contain
## this subplot
matchidx <- as.data.frame( which( layout == ii,
arr.ind = TRUE ) )
print( plots[[ ii ]], vp = viewport(
layout.pos.row = matchidx$row,
layout.pos.col = matchidx$col ) )
}
}
}
invisible( last_plot() )
}
##' @title Plot xts class data
##' @description Plotting a \pkg{xts} time series with \pkg{ggplot2}
##' in a convenient format.
##'
##' @details Plots all objects of class \pkg{xts}. With the main
##' argument a title can be provided.
##'
##' @param data.input Class \pkg{xts} time series, which should be
##' visualized.
##' @param main Title of the plot to generate. Default: "Time series"
##' @param ylab Label of the y-axis of the plot to generate. Default =
##' NULL
##' @param ... Additional parameters for the multiplot function.
##'
##' @family plot
##' @export
##' @importFrom xts is.xts
##' @import ggplot2
##' @return Nothing. The multiplot function is run in the last step.
##' @author Philipp Mueller
ttplot <- function( data.input, main = "Time series", ylab = NULL ){
if ( !is.xts( data.input ) )
stop( "Please supply an object of class 'xts'!" )
if ( is.null( ylab ) )
ylab <- names( data.input )
if ( is.null( ylab ) )
ylab <- ""
## Checks if the time series has daily values or if it is blocked.
## In the latter case the create values can not be provided as the
## input of the x axis because ggplot2 does not know how to scale
## them.
x.time.unit <- index( data.input )[[ 2 ]] -
index( data.input )[[ 1 ]]
x.df <- data.frame( date = as.numeric( index( data.input ) ),
value = as.numeric( data.input ) )
ggplot( data = x.df, aes( x = date, y = value ) ) +
geom_point( colour = "darkorange" ) +
geom_line( colour = "navy" ) +
ggtitle( main ) + xlab( "Time" ) + ylab( ylab ) +
theme_bw() +
theme( panel.grid.major = element_line( colour = "grey75" ),
panel.grid.minor = element_line( colour = "grey80" ) )
return( last_plot() )
}
##' @title Plots the GEV likelihood
##' @description Generates three 2D slices of the likelihood of a
##' provided time series.
##'
##' @details The likelihood will be plotted around its optimal
##' value. To determine it the time series will be fitted via
##' \code{\link{fit.gev}}. Its also possible to provide another
##' point the likelihood will be centered around using
##' \strong{center}.
##'
##' The likelihood will be calculated using the C++-based function
##' \strong{likelihood_GEV}.
##'
##' @param time.series Data, for which the likelihood function is
##' going to be calculated.
##' @param location.lim Boundaries of the plotted likelihood
##' region. Default = +/- 2.5 times the minimum value.
##' @param scale.lim Boundaries of the plotted likelihood
##' region. Default = +/- 2.5 times the minimum value.
##' @param shape.lim Boundaries of the plotted likelihood
##' region. Default = +/- 2.5 times the minimum value.
##' @param center It set this point (3D) will be the new center of the
##' likelihood plot. Default = NULL
##' @param initial Initial parameters which can be provided for the
##' fit of the time series.
##' @param main Title for the multiplot. Default = NULL.
##' @param true.minima Provide the true minima of the neg
##' log-likelihood function if \code{\link{fit.gev}} does not
##' provide it in the first run. Default = NULL.
##'
##' @import ggplot2
##' @return returns a multiplot of planes through the negative
##' log-likelihood of the GEV function calculated using the input
##' time series.
##' @author Philipp Mueller
likelihood.plot <- function( time.series, location.lim = NULL,
scale.lim = NULL, shape.lim = NULL,
center = NULL, initial = NULL,
main = NULL, true.minima = NULL ){
## determining the center and the limits
if ( is.null( true.minima ) ){
time.series.mle <- fit.gev( x = time.series,
initial = initial,
silent = TRUE )$par
} else
time.series.mle <- true.minima
if ( is.null( center ) )
center <- time.series.mle
if ( is.null( location.lim ) )
location.lim <- c( center[ 1 ] - time.series.mle[ 1 ]* .5,
center[ 1 ] + time.series.mle[ 1 ]* .5 )
if ( is.null( scale.lim ) )
scale.lim <- c( center[ 2 ] - time.series.mle[ 2 ]* .5,
center[ 2 ] + time.series.mle[ 2 ]* .5 )
if( scale.lim[ 1 ] < 0 )
scale.lim[ 1 ] <- 0
if ( is.null( shape.lim ) )
shape.lim <- c( center[ 3 ] - time.series.mle[ 3 ]* .5,
center[ 3 ] + time.series.mle[ 3 ]* .5 )
## 2D likelihood
## some more point because its just the plane this time
number.of.points <- 200
## setting the limits of the scale and shape parameter if not
## provided yet
location.range <- seq( location.lim[ 1 ], location.lim[ 2 ], ,
number.of.points )
scale.range <- seq( scale.lim[ 1 ], scale.lim[ 2 ], ,
number.of.points )
shape.range <- seq( shape.lim[ 1 ], shape.lim[ 2 ], ,
number.of.points )
loc.sc.plane <- data.frame(
expand.grid( location.range, scale.range )[ , 1 ],
expand.grid( location.range, scale.range )[ , 2 ],
rep( center[ 3 ], number.of.points^ 2 ) )
loc.sh.plane <- data.frame(
expand.grid( location.range, shape.range )[ , 1 ],
rep( center[ 2 ], number.of.points^ 2 ),
expand.grid( location.range, shape.range )[ , 2 ] )
sc.sh.plane <- data.frame(
rep( center[ 1 ], number.of.points^ 2 ),
expand.grid( scale.range, shape.range )[ , 1 ],
expand.grid( scale.range, shape.range )[ , 2 ] )
plot.plane <- function( par.plane, par.plane.names,
time.series = time.series ){
names( par.plane ) <- c( "location", "scale", "shape" )
likelihood.plane <- .Call( 'likelihood_GEV', PACKAGE = 'climex',
par.plane, time.series )
## Maybe introducing a cutoff in the likelihood values make the
## contour plot work again
likelihood.plane[ likelihood.plane > 1E8 ] <- 1E8
min.likelihood <- min( likelihood.plane, na.rm = TRUE )
if ( min.likelihood < 0 ){
max.likelihood.low <- -min.likelihood* 1.5
## comment the next line to have a cutoff at 0 in the likelihood
## surface plot
likelihood.plane <- likelihood.plane - min.likelihood + 0.01
warning( "In some regions the likelihood is negative. An offset is used in order to still plot it in a logarithmic way" )
} else
max.likelihood.low <- min.likelihood* 1.5
par.plane$likelihood <- par.plane$likelihood.low <-
likelihood.plane
par.plane$likelihood.low[ par.plane$likelihood.low >
max.likelihood.low ] <-
max.likelihood.low
names( par.plane ) <- c( par.plane.names, "likelihood.low",
"likelihood" )
gg.plane <-
ggplot() + geom_raster( data = par.plane,
aes( x = x, y = y, fill = likelihood ),
na.rm = TRUE ) +
geom_contour( data = par.plane,
aes( x = x, y = y, z = likelihood.low ),
colour = grDevices::rgb( 1, .55, 0 ),
na.rm = TRUE ) +
scale_fill_gradientn(
colours = rev( RColorBrewer::brewer.pal( 7, "Blues" ) ),
na.value = "white", trans = "log" ) +
theme_bw()
return( last_plot() )
}
plot.center <- data.frame( x = center[ 1 ], y = center[ 2 ],
z = center[ 3 ] )
multiplot( plotlist = list(
plot.plane( loc.sh.plane,
par.plane.names = c( "x", "a", "y" ),
time.series ) +
geom_point( data = plot.center, aes( x = x, y = z ),
shape = 4 ) +
xlab( "location" ) + ylab( "shape" ),
plot.plane( sc.sh.plane,
par.plane.names = c( "a", "x", "y" ),
time.series ) +
geom_point( data = plot.center, aes( x = y, y = z ),
shape = 4 ) +
xlab( "scale" ) + ylab( "shape" ),
plot.plane( loc.sc.plane,
par.plane.names = c( "x", "y", "a" ),
time.series ) +
geom_point( data = plot.center, aes( x = x, y = y ),
shape = 4 ) +
xlab( "location" ) + ylab( "scale" ) ),
main = main )
}
##' @title Plots a GEV fit
##' @description Plots the GEV function fitted using
##' \code{\link{fit.gev}}
##'
##' @details Uses \pkg{ggplot2}.
##'
##' Since this function is also uses in the shiny app
##' \emph{climexUI}, the number of displayed bins has to be
##' adjusted.
##'
##' @param x Fitted GEV object of class \emph{climex.fit.gev}.
##' @param bin.factor Multiplying the length of \strong{x} by this
##' factor, results in the number of bins used in this plot. Default
##' = NULL
##' @param ... Additional parameters. They won't be handled in the
##' function. This argument is only present to ensure S3 generic
##' consistency with respect to the \code{\link[graphics]{plot}}
##' function.
##'
##' @export
##' @import ggplot2
##' @return ggplot2 object.
plot.climex.fit.gev <- function( x, bin.factor = NULL, ... ){
x.data <- x$x
if ( is.null( bin.factor ) ){
bin.factor <- ( ( ( length( x.data ) - 1 )*100/
length( x.data ) ) %/% 5 )* 0.025
}
bin.width <- ( ( max( x.data ) - min( x.data ) )/
( bin.factor * length( x.data ) ) )* 1.1
## Range of the supplied data
x.lim <- c( min( x.data, na.rm = TRUE ),
max( x.data, na.rm = TRUE ) )
## The the calculated density falls below this threshold it won't
## get plotted anymore.
threshold.pdf.plot <- 1E-4
if ( x$control$extreme.type == "min" ){
x.lim.rev <- -1 * rev( x.lim )
## In the beginning I used the coefficients of the GEV distribution
## to determine the range of evaluation of the density function. But
## this will fail for distributions with bigger absolute shape
## value. Instead I will use the range of the supplied data
plot.gev.range <- seq( x.lim.rev[ 1 ] - x$par[ 2 ]* 10,
x.lim.rev[ 2 ] + x$par[ 2 ]* 10, 0.01 )
plot.gev.data <- data.frame(
x.plot = -1* plot.gev.range,
y.plot = gev.density( x$par * c( -1, 1, 1 ),
plot.gev.range ) )
} else if ( x$control$extreme.type == "max" ){
## In the beginning I used the coefficients of the GEV distribution
## to determine the range of evaluation of the density function. But
## this will fail for distributions with bigger absolute shape
## value. Instead I will use the range of the supplied data
plot.gev.range <- seq( x.lim[ 1 ] - x$par[ 2 ]* 10,
x.lim[ 2 ] + x$par[ 2 ]* 10, 0.01 )
plot.gev.data <- data.frame(
x.plot = plot.gev.range,
y.plot = gev.density( x$par, plot.gev.range ) )
} else {
stop( "Unexpected extreme.type in the provided fitting object" )
}
## Removing all NaN to have less points and warnings
plot.gev.data <- plot.gev.data[
!is.nan( plot.gev.data$y.plot ), ]
## Cutting of the density at values below the threshold in the tails
## Lower tail
plot.gev.data <- plot.gev.data[
( plot.gev.data$x.plot >
plot.gev.data$x.plot[ which.max( plot.gev.data$y.plot ) ]
) | ( plot.gev.data$y.plot > threshold.pdf.plot ), ]
## Upper tail
plot.gev.data <- plot.gev.data[
( plot.gev.data$x.plot <
plot.gev.data$x.plot[ which.max( plot.gev.data$y.plot ) ]
) | ( plot.gev.data$y.plot > threshold.pdf.plot ), ]
## If the first and last entry is not a zero, add one!
## Else the density won't look like a density at all.
if ( plot.gev.data$y.plot[ 1 ] > 0 ){
plot.gev.data <- data.frame(
x.plot = c( 2* plot.gev.data$x.plot[ 1 ] -
plot.gev.data$x.plot[ 2 ],
plot.gev.data$x.plot ),
y.plot = c( 0, plot.gev.data$y.plot ) )
}
if ( plot.gev.data$y.plot[ nrow( plot.gev.data ) ] > 0 ){
plot.gev.data[ nrow( plot.gev.data ) + 1, ] <-
c( 2* plot.gev.data$x.plot[ nrow( plot.gev.data ) ] -
plot.gev.data$x.plot[ nrow( plot.gev.data ) - 1 ],
0 )
}
plot.gev.lim <- c( min( plot.gev.data$x.plot ),
max( plot.gev.data$x.plot ) )
## Determining the limits of the actual plot
if ( plot.gev.lim[ 1 ] > x.lim[ 1 ] )
plot.gev.lim[ 1 ] <- x.lim[ 1 ] - abs( x.lim[ 1 ] )* 0.05
if ( plot.gev.lim[ 2 ] < x.lim[ 2 ] )
plot.gev.lim[ 2 ] <- x.lim[ 2 ] + abs( x.lim[ 2 ] )* 0.05
plot.lim <- c( min( plot.gev.lim[ 1 ],
min( x.data ) - bin.width ),
max( plot.gev.lim[ 2 ],
max( x.data) + bin.width ) )
## Whenever the density gets bigger than one, cut it
if ( max( plot.gev.data$y.plot ) > 1 ){
y.lim <- c( 0, ( max( stats::density( x.data )$y )* 1.5 ) )
} else {
y.lim <- NULL
}
ggplot() + geom_histogram(
data = data.frame( x = x.data ),
colour = grDevices::rgb( .098, .098, .44 ),
alpha = 1,
bins = bin.factor* length( x.data ),
aes( x = x, y = ..density..,
fill = "#7171EC" ) ) +
geom_polygon( data = plot.gev.data, alpha = 0.7,
colour = grDevices::rgb( .098, .098, .44 ),
aes( x = x.plot, y = y.plot,
fill = grDevices::rgb( 1, .55, 0 ) ) ) +
scale_fill_manual( values = c( "#7171EC",
grDevices::rgb( 1, .55, 0 ) ),
labels = c( "Histogram", "Fitted GEV" ) ) +
theme_bw() + ylab( "Density" ) +
coord_cartesian( xlim = plot.lim, ylim = y.lim ) +
theme( legend.title = element_blank() )
return( last_plot() )
}
##' @title Plot a GP fit
##' @description Plots the GP function fitted using
##' \code{\link{fit.gpd}}
##'
##' @details Uses \pkg{ggplot2}.
##'
##' Since this function is also uses in the shiny app
##' \emph{climexUI}, the number of displayed bins has to be
##' adjusted.
##'
##' @param x Fitted GP object of class \emph{climex.fit.gpd}.
##' @param bin.factor Multiplying the length of \strong{x} by this
##' factor, results in the number of bins used in this plot. Default
##' = NULL
##' @param ... Additional parameters. They won't be handled in the
##' function. This argument is only present to ensure S3 generic
##' consistency with respect to the \code{\link[graphics]{plot}}
##' function.
##'
##' @export
##' @import ggplot2
##' @return ggplot2 object.
plot.climex.fit.gpd <- function( x, bin.factor = NULL, ... ){
if ( is.null( x$threshold ) ){
stop( "Please provide a threshold argument to fit.gpd() in order to plot the result" )
}
x.data <- as.numeric( x$x + x$threshold )
if ( is.null( bin.factor ) ){
bin.factor <- ( ( ( length( x.data ) - 1 )*100/
length( x.data ) ) %/% 5 )* 0.025
}
bin.width <- ( ( max( x.data ) - min( x.data ) )/
( bin.factor * length( x.data ) ) )* 1.1
## Range of the supplied data
x.lim <- c( min( x.data, na.rm = TRUE ),
max( x.data, na.rm = TRUE ) )
## The the calculated density falls below this threshold it won't
## get plotted anymore.
threshold.pdf.plot <- 1E-4
if ( x$control$extreme.type == "min" ){
x.lim.rev <- -1 * rev( x.lim )
## In the beginning I used the coefficients of the GP distribution
## to determine the range of evaluation of the density function. But
## this will fail for distributions with bigger absolute shape
## value. Instead I will use the range of the supplied data
plot.range <- seq( x.lim.rev[ 1 ] - bin.width,
x.lim.rev[ 2 ] + x$par[ 1 ]* 10, 0.01 )
plot.data <- data.frame(
x.plot = -1* plot.range,
y.plot = gpd.density( x$par, x$threshold * -1,
plot.range ) )
} else if ( x$control$extreme.type == "max" ){
## In the beginning I used the coefficients of the GP distribution
## to determine the range of evaluation of the density function. But
## this will fail for distributions with bigger absolute shape
## value. Instead I will use the range of the supplied data
plot.range <- seq( x.lim[ 1 ] - bin.width,
x.lim[ 2 ] + x$par[ 1 ]* 10, 0.01 )
plot.data <- data.frame(
x.plot = plot.range,
y.plot = gpd.density( x$par, x$threshold,
plot.range ) )
} else {
stop( "Unexpected extreme.type in provided object" )
}
## Removing all NaN to have less points and warnings
plot.data <- plot.data[
!is.nan( plot.data$y.plot ), ]
## Cutting of the density at values below the threshold in the tails
## Lower tail
plot.data <- plot.data[
( plot.data$x.plot >
plot.data$x.plot[ which.max( plot.data$y.plot ) ]
) | ( plot.data$y.plot > threshold.pdf.plot ), ]
## Upper tail
plot.data <- plot.data[
( plot.data$x.plot <
plot.data$x.plot[ which.max( plot.data$y.plot ) ]
) | ( plot.data$y.plot > threshold.pdf.plot ), ]
## If the first and last entry is not a zero, add one!
## Else the density won't look like a density at all.
if ( plot.data$y.plot[ 1 ] > 0 ){
plot.data <- data.frame(
x.plot = c( 2* plot.data$x.plot[ 1 ] -
plot.data$x.plot[ 2 ],
plot.data$x.plot ),
y.plot = c( 0, plot.data$y.plot ) )
}
if ( plot.data$y.plot[ nrow( plot.data ) ] > 0 ){
plot.data[ nrow( plot.data ) + 1, ] <-
c( 2* plot.data$x.plot[ nrow( plot.data ) ] -
plot.data$x.plot[ nrow( plot.data ) - 1 ],
0 )
}
plot.lim <- c( min( plot.data$x.plot ),
max( plot.data$x.plot ) )
## Using the threshold value defined above the range of
## the density plot will be determined
plot.lim <- c( plot.data[[ 1 ]][
which.min( abs( plot.data[[ 2 ]][
1 : which.max( plot.data[[ 2 ]] ) ] -
threshold.pdf.plot ) ) ],
plot.data[[ 1 ]][ which.min( (
plot.data[[ 2 ]][ which.max( plot.data[[ 2 ]] )
: length( plot.data[[ 2 ]] ) ] -
threshold.pdf.plot ) ) +
which.max( plot.data[[ 2 ]] ) - 1 ] )
## Determining the limits of the actual plot
if ( plot.lim[ 1 ] > x.lim[ 1 ] )
plot.lim[ 1 ] <- x.lim[ 1 ] - abs( x.lim[ 1 ] )* 0.05
if ( plot.lim[ 2 ] < x.lim[ 2 ] )
plot.lim[ 2 ] <- x.lim[ 2 ] + abs( x.lim[ 2 ] )* 0.05
plot.lim <- c( min( plot.lim[ 1 ],
min( x.data ) - bin.width ),
max( plot.lim[ 2 ],
max( x.data) + bin.width ) )
## Whenever the density gets bigger than one, cut it
if ( max( plot.data$y.plot ) > 1 ){
y.lim <- c( 0, ( max( stats::density( x.data )$y )* 1.5 ) )
} else {
y.lim <- NULL
}
ggplot() + geom_histogram(
data = data.frame( x = x.data ),
colour = grDevices::rgb( .098, .098, .44 ),
alpha = 1,
bins = bin.factor* length( x.data ),
aes( x = x.data, y = ..density..,
fill = "#7171EC" ) ) +
geom_polygon( data = plot.data, alpha = 0.7,
colour = grDevices::rgb( .098, .098, .44 ),
aes( x = x.plot, y = y.plot,
fill = grDevices::rgb( 1, .55, 0 ) ) ) +
scale_fill_manual( values = c( "#7171EC",
grDevices::rgb( 1, .55, 0 ) ),
labels = c( "Histogram", "Fitted GP" ) ) +
theme_bw() + ylab( "Density" ) +
coord_cartesian( xlim = plot.lim, ylim = y.lim ) +
theme( legend.title = element_blank() )
return( last_plot() )
}
## End of plot.R
philmu/climex documentation built on July 11, 2022, 3:23 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions plot.climex.fit.gpd plot.climex.fit.gev likelihood.plot ttplot multiplot
Documented in likelihood.plot multiplot plot.climex.fit.gev plot.climex.fit.gpd ttplot
### plot.R - Functions for plotting either the 'xts' class data or the
## fitting results
##' @title Plot multiple time series.
##' @description This function takes several \pkg{ggplot2} plots and
##' arranges them on one page.
##' @details Well, it is not really related to extreme value fitting
##' and I am aware of the fact that its quite bad style to export
##' auxiliary functions. But this one is just so handy. If the
##' \strong{layout} is something like
##' \code{matrix( c( 1, 2, 3, 3 ), nrow = 2, byrow = TRUE )},
##' then plot 1 will go in the upper left, 2 will go in the upper
##' right, and 3 will go all the way across the bottom. ggplot
##' objects can be passed in ..., or to \strong{plotlist} (as a list
## of \emph{ggplot} objects).
##'
##' @param plotlist Alternative way to pass \emph{ggplot} objects.
##' @param tt.title Same as \strong{main}.
##' @param main Title of the overall plot.
##' @param ... Additional arguments for ggplot.
##' @param cols Number of columns in layout. Default = 1.
##' @param layout A matrix specifying the layout. If present,
##' \strong{cols} is ignored.
##'
##' @family plot
##' @author Paul Teetor (original), Philpp Mueller (minor changes)
##' @examples
##' anomalies.potsdam <- anomalies( temp.potsdam )
##' multiplot( anomalies.potsdam, temp.potsdam, cols = 2, main =
##' "Difference between the pure data of the daily maximum
##' temperatures of the Potsdam station and their anomalies" )
##' @return print object.
##' @export
##' @import ggplot2
##' @import grid
multiplot <- function( plotlist = NULL, tt.title = main, main = NULL,
..., cols = 1, layout = NULL ) {
## Make a list from the ... arguments and plotlist
plots <- c( list( ... ), plotlist )
numPlots = length( plots )
## If layout is NULL, then use 'cols' to determine layout
if ( is.null( layout ) ) {
## Make the panel
## ncol: Number of columns of plots
## nrow: Number of rows needed, calculated from # of cols
layout <- matrix( seq( 1, cols* ceiling ( numPlots/ cols ) ),
ncol = cols, nrow = ceiling( numPlots/ cols ) )
}
if ( numPlots == 1 ) {
print( plots[[ 1 ]] )
} else {
## Set up the page
grid.newpage()
if ( !is.null( tt.title ) ){
pushViewport( viewport(
layout = grid.layout(
nrow( layout ) + 1, ncol( layout ),
heights = unit( rep( 1, nrow( layout ) + 1 ),
c( "lines", rep( "null",
nrow( layout) ) ) ) ) ) )
} else
pushViewport( viewport(
layout = grid.layout( nrow( layout ), ncol( layout ) ) ) )
## Make each plot, in the correct location
if ( !is.null( tt.title ) ){
grid.text( tt.title, vp = viewport(
layout.pos.row = 1,
layout.pos.col = c( 1, ncol( layout ) )
) )
for ( ii in 1 : numPlots ) {
## Get the i,j matrix positions of the regions that contain
## this subplot
matchidx <- as.data.frame( which( layout == ii,
arr.ind = TRUE ) )
print( plots[[ ii ]], vp = viewport(
layout.pos.row = matchidx$row + 1,
layout.pos.col = matchidx$col ) )
}
} else {
for ( ii in 1 : numPlots ) {
## Get the i,j matrix positions of the regions that contain
## this subplot
matchidx <- as.data.frame( which( layout == ii,
arr.ind = TRUE ) )
print( plots[[ ii ]], vp = viewport(
layout.pos.row = matchidx$row,
layout.pos.col = matchidx$col ) )
}
}
}
invisible( last_plot() )
}
##' @title Plot xts class data
##' @description Plotting a \pkg{xts} time series with \pkg{ggplot2}
##' in a convenient format.
##'
##' @details Plots all objects of class \pkg{xts}. With the main
##' argument a title can be provided.
##'
##' @param data.input Class \pkg{xts} time series, which should be
##' visualized.
##' @param main Title of the plot to generate. Default: "Time series"
##' @param ylab Label of the y-axis of the plot to generate. Default =
##' NULL
##' @param ... Additional parameters for the multiplot function.
##'
##' @family plot
##' @export
##' @importFrom xts is.xts
##' @import ggplot2
##' @return Nothing. The multiplot function is run in the last step.
##' @author Philipp Mueller
ttplot <- function( data.input, main = "Time series", ylab = NULL ){
if ( !is.xts( data.input ) )
stop( "Please supply an object of class 'xts'!" )
if ( is.null( ylab ) )
ylab <- names( data.input )
if ( is.null( ylab ) )
ylab <- ""
## Checks if the time series has daily values or if it is blocked.
## In the latter case the create values can not be provided as the
## input of the x axis because ggplot2 does not know how to scale
## them.
x.time.unit <- index( data.input )[[ 2 ]] -
index( data.input )[[ 1 ]]
x.df <- data.frame( date = as.numeric( index( data.input ) ),
value = as.numeric( data.input ) )
ggplot( data = x.df, aes( x = date, y = value ) ) +
geom_point( colour = "darkorange" ) +
geom_line( colour = "navy" ) +
ggtitle( main ) + xlab( "Time" ) + ylab( ylab ) +
theme_bw() +
theme( panel.grid.major = element_line( colour = "grey75" ),
panel.grid.minor = element_line( colour = "grey80" ) )
return( last_plot() )
}
##' @title Plots the GEV likelihood
##' @description Generates three 2D slices of the likelihood of a
##' provided time series.
##'
##' @details The likelihood will be plotted around its optimal
##' value. To determine it the time series will be fitted via
##' \code{\link{fit.gev}}. Its also possible to provide another
##' point the likelihood will be centered around using
##' \strong{center}.
##'
##' The likelihood will be calculated using the C++-based function
##' \strong{likelihood_GEV}.
##'
##' @param time.series Data, for which the likelihood function is
##' going to be calculated.
##' @param location.lim Boundaries of the plotted likelihood
##' region. Default = +/- 2.5 times the minimum value.
##' @param scale.lim Boundaries of the plotted likelihood
##' region. Default = +/- 2.5 times the minimum value.
##' @param shape.lim Boundaries of the plotted likelihood
##' region. Default = +/- 2.5 times the minimum value.
##' @param center It set this point (3D) will be the new center of the
##' likelihood plot. Default = NULL
##' @param initial Initial parameters which can be provided for the
##' fit of the time series.
##' @param main Title for the multiplot. Default = NULL.
##' @param true.minima Provide the true minima of the neg
##' log-likelihood function if \code{\link{fit.gev}} does not
##' provide it in the first run. Default = NULL.
##'
##' @import ggplot2
##' @return returns a multiplot of planes through the negative
##' log-likelihood of the GEV function calculated using the input
##' time series.
##' @author Philipp Mueller
likelihood.plot <- function( time.series, location.lim = NULL,
scale.lim = NULL, shape.lim = NULL,
center = NULL, initial = NULL,
main = NULL, true.minima = NULL ){
## determining the center and the limits
if ( is.null( true.minima ) ){
time.series.mle <- fit.gev( x = time.series,
initial = initial,
silent = TRUE )$par
} else
time.series.mle <- true.minima
if ( is.null( center ) )
center <- time.series.mle
if ( is.null( location.lim ) )
location.lim <- c( center[ 1 ] - time.series.mle[ 1 ]* .5,
center[ 1 ] + time.series.mle[ 1 ]* .5 )
if ( is.null( scale.lim ) )
scale.lim <- c( center[ 2 ] - time.series.mle[ 2 ]* .5,
center[ 2 ] + time.series.mle[ 2 ]* .5 )
if( scale.lim[ 1 ] < 0 )
scale.lim[ 1 ] <- 0
if ( is.null( shape.lim ) )
shape.lim <- c( center[ 3 ] - time.series.mle[ 3 ]* .5,
center[ 3 ] + time.series.mle[ 3 ]* .5 )
## 2D likelihood
## some more point because its just the plane this time
number.of.points <- 200
## setting the limits of the scale and shape parameter if not
## provided yet
location.range <- seq( location.lim[ 1 ], location.lim[ 2 ], ,
number.of.points )
scale.range <- seq( scale.lim[ 1 ], scale.lim[ 2 ], ,
number.of.points )
shape.range <- seq( shape.lim[ 1 ], shape.lim[ 2 ], ,
number.of.points )
loc.sc.plane <- data.frame(
expand.grid( location.range, scale.range )[ , 1 ],
expand.grid( location.range, scale.range )[ , 2 ],
rep( center[ 3 ], number.of.points^ 2 ) )
loc.sh.plane <- data.frame(
expand.grid( location.range, shape.range )[ , 1 ],
rep( center[ 2 ], number.of.points^ 2 ),
expand.grid( location.range, shape.range )[ , 2 ] )
sc.sh.plane <- data.frame(
rep( center[ 1 ], number.of.points^ 2 ),
expand.grid( scale.range, shape.range )[ , 1 ],
expand.grid( scale.range, shape.range )[ , 2 ] )
plot.plane <- function( par.plane, par.plane.names,
time.series = time.series ){
names( par.plane ) <- c( "location", "scale", "shape" )
likelihood.plane <- .Call( 'likelihood_GEV', PACKAGE = 'climex',
par.plane, time.series )
## Maybe introducing a cutoff in the likelihood values make the
## contour plot work again
likelihood.plane[ likelihood.plane > 1E8 ] <- 1E8
min.likelihood <- min( likelihood.plane, na.rm = TRUE )
if ( min.likelihood < 0 ){
max.likelihood.low <- -min.likelihood* 1.5
## comment the next line to have a cutoff at 0 in the likelihood
## surface plot
likelihood.plane <- likelihood.plane - min.likelihood + 0.01
warning( "In some regions the likelihood is negative. An offset is used in order to still plot it in a logarithmic way" )
} else
max.likelihood.low <- min.likelihood* 1.5
par.plane$likelihood <- par.plane$likelihood.low <-
likelihood.plane
par.plane$likelihood.low[ par.plane$likelihood.low >
max.likelihood.low ] <-
max.likelihood.low
names( par.plane ) <- c( par.plane.names, "likelihood.low",
"likelihood" )
gg.plane <-
ggplot() + geom_raster( data = par.plane,
aes( x = x, y = y, fill = likelihood ),
na.rm = TRUE ) +
geom_contour( data = par.plane,
aes( x = x, y = y, z = likelihood.low ),
colour = grDevices::rgb( 1, .55, 0 ),
na.rm = TRUE ) +
scale_fill_gradientn(
colours = rev( RColorBrewer::brewer.pal( 7, "Blues" ) ),
na.value = "white", trans = "log" ) +
theme_bw()
return( last_plot() )
}
plot.center <- data.frame( x = center[ 1 ], y = center[ 2 ],
z = center[ 3 ] )
multiplot( plotlist = list(
plot.plane( loc.sh.plane,
par.plane.names = c( "x", "a", "y" ),
time.series ) +
geom_point( data = plot.center, aes( x = x, y = z ),
shape = 4 ) +
xlab( "location" ) + ylab( "shape" ),
plot.plane( sc.sh.plane,
par.plane.names = c( "a", "x", "y" ),
time.series ) +
geom_point( data = plot.center, aes( x = y, y = z ),
shape = 4 ) +
xlab( "scale" ) + ylab( "shape" ),
plot.plane( loc.sc.plane,
par.plane.names = c( "x", "y", "a" ),
time.series ) +
geom_point( data = plot.center, aes( x = x, y = y ),
shape = 4 ) +
xlab( "location" ) + ylab( "scale" ) ),
main = main )
}
##' @title Plots a GEV fit
##' @description Plots the GEV function fitted using
##' \code{\link{fit.gev}}
##'
##' @details Uses \pkg{ggplot2}.
##'
##' Since this function is also uses in the shiny app
##' \emph{climexUI}, the number of displayed bins has to be
##' adjusted.
##'
##' @param x Fitted GEV object of class \emph{climex.fit.gev}.
##' @param bin.factor Multiplying the length of \strong{x} by this
##' factor, results in the number of bins used in this plot. Default
##' = NULL
##' @param ... Additional parameters. They won't be handled in the
##' function. This argument is only present to ensure S3 generic
##' consistency with respect to the \code{\link[graphics]{plot}}
##' function.
##'
##' @export
##' @import ggplot2
##' @return ggplot2 object.
plot.climex.fit.gev <- function( x, bin.factor = NULL, ... ){
x.data <- x$x
if ( is.null( bin.factor ) ){
bin.factor <- ( ( ( length( x.data ) - 1 )*100/
length( x.data ) ) %/% 5 )* 0.025
}
bin.width <- ( ( max( x.data ) - min( x.data ) )/
( bin.factor * length( x.data ) ) )* 1.1
## Range of the supplied data
x.lim <- c( min( x.data, na.rm = TRUE ),
max( x.data, na.rm = TRUE ) )
## The the calculated density falls below this threshold it won't
## get plotted anymore.
threshold.pdf.plot <- 1E-4
if ( x$control$extreme.type == "min" ){
x.lim.rev <- -1 * rev( x.lim )
## In the beginning I used the coefficients of the GEV distribution
## to determine the range of evaluation of the density function. But
## this will fail for distributions with bigger absolute shape
## value. Instead I will use the range of the supplied data
plot.gev.range <- seq( x.lim.rev[ 1 ] - x$par[ 2 ]* 10,
x.lim.rev[ 2 ] + x$par[ 2 ]* 10, 0.01 )
plot.gev.data <- data.frame(
x.plot = -1* plot.gev.range,
y.plot = gev.density( x$par * c( -1, 1, 1 ),
plot.gev.range ) )
} else if ( x$control$extreme.type == "max" ){
## In the beginning I used the coefficients of the GEV distribution
## to determine the range of evaluation of the density function. But
## this will fail for distributions with bigger absolute shape
## value. Instead I will use the range of the supplied data
plot.gev.range <- seq( x.lim[ 1 ] - x$par[ 2 ]* 10,
x.lim[ 2 ] + x$par[ 2 ]* 10, 0.01 )
plot.gev.data <- data.frame(
x.plot = plot.gev.range,
y.plot = gev.density( x$par, plot.gev.range ) )
} else {
stop( "Unexpected extreme.type in the provided fitting object" )
}
## Removing all NaN to have less points and warnings
plot.gev.data <- plot.gev.data[
!is.nan( plot.gev.data$y.plot ), ]
## Cutting of the density at values below the threshold in the tails
## Lower tail
plot.gev.data <- plot.gev.data[
( plot.gev.data$x.plot >
plot.gev.data$x.plot[ which.max( plot.gev.data$y.plot ) ]
) | ( plot.gev.data$y.plot > threshold.pdf.plot ), ]
## Upper tail
plot.gev.data <- plot.gev.data[
( plot.gev.data$x.plot <
plot.gev.data$x.plot[ which.max( plot.gev.data$y.plot ) ]
) | ( plot.gev.data$y.plot > threshold.pdf.plot ), ]
## If the first and last entry is not a zero, add one!
## Else the density won't look like a density at all.
if ( plot.gev.data$y.plot[ 1 ] > 0 ){
plot.gev.data <- data.frame(
x.plot = c( 2* plot.gev.data$x.plot[ 1 ] -
plot.gev.data$x.plot[ 2 ],
plot.gev.data$x.plot ),
y.plot = c( 0, plot.gev.data$y.plot ) )
}
if ( plot.gev.data$y.plot[ nrow( plot.gev.data ) ] > 0 ){
plot.gev.data[ nrow( plot.gev.data ) + 1, ] <-
c( 2* plot.gev.data$x.plot[ nrow( plot.gev.data ) ] -
plot.gev.data$x.plot[ nrow( plot.gev.data ) - 1 ],
0 )
}
plot.gev.lim <- c( min( plot.gev.data$x.plot ),
max( plot.gev.data$x.plot ) )
## Determining the limits of the actual plot
if ( plot.gev.lim[ 1 ] > x.lim[ 1 ] )
plot.gev.lim[ 1 ] <- x.lim[ 1 ] - abs( x.lim[ 1 ] )* 0.05
if ( plot.gev.lim[ 2 ] < x.lim[ 2 ] )
plot.gev.lim[ 2 ] <- x.lim[ 2 ] + abs( x.lim[ 2 ] )* 0.05
plot.lim <- c( min( plot.gev.lim[ 1 ],
min( x.data ) - bin.width ),
max( plot.gev.lim[ 2 ],
max( x.data) + bin.width ) )
## Whenever the density gets bigger than one, cut it
if ( max( plot.gev.data$y.plot ) > 1 ){
y.lim <- c( 0, ( max( stats::density( x.data )$y )* 1.5 ) )
} else {
y.lim <- NULL
}
ggplot() + geom_histogram(
data = data.frame( x = x.data ),
colour = grDevices::rgb( .098, .098, .44 ),
alpha = 1,
bins = bin.factor* length( x.data ),
aes( x = x, y = ..density..,
fill = "#7171EC" ) ) +
geom_polygon( data = plot.gev.data, alpha = 0.7,
colour = grDevices::rgb( .098, .098, .44 ),
aes( x = x.plot, y = y.plot,
fill = grDevices::rgb( 1, .55, 0 ) ) ) +
scale_fill_manual( values = c( "#7171EC",
grDevices::rgb( 1, .55, 0 ) ),
labels = c( "Histogram", "Fitted GEV" ) ) +
theme_bw() + ylab( "Density" ) +
coord_cartesian( xlim = plot.lim, ylim = y.lim ) +
theme( legend.title = element_blank() )
return( last_plot() )
}
##' @title Plot a GP fit
##' @description Plots the GP function fitted using
##' \code{\link{fit.gpd}}
##'
##' @details Uses \pkg{ggplot2}.
##'
##' Since this function is also uses in the shiny app
##' \emph{climexUI}, the number of displayed bins has to be
##' adjusted.
##'
##' @param x Fitted GP object of class \emph{climex.fit.gpd}.
##' @param bin.factor Multiplying the length of \strong{x} by this
##' factor, results in the number of bins used in this plot. Default
##' = NULL
##' @param ... Additional parameters. They won't be handled in the
##' function. This argument is only present to ensure S3 generic
##' consistency with respect to the \code{\link[graphics]{plot}}
##' function.
##'
##' @export
##' @import ggplot2
##' @return ggplot2 object.
plot.climex.fit.gpd <- function( x, bin.factor = NULL, ... ){
if ( is.null( x$threshold ) ){
stop( "Please provide a threshold argument to fit.gpd() in order to plot the result" )
}
x.data <- as.numeric( x$x + x$threshold )
if ( is.null( bin.factor ) ){
bin.factor <- ( ( ( length( x.data ) - 1 )*100/
length( x.data ) ) %/% 5 )* 0.025
}
bin.width <- ( ( max( x.data ) - min( x.data ) )/
( bin.factor * length( x.data ) ) )* 1.1
## Range of the supplied data
x.lim <- c( min( x.data, na.rm = TRUE ),
max( x.data, na.rm = TRUE ) )
## The the calculated density falls below this threshold it won't
## get plotted anymore.
threshold.pdf.plot <- 1E-4
if ( x$control$extreme.type == "min" ){
x.lim.rev <- -1 * rev( x.lim )
## In the beginning I used the coefficients of the GP distribution
## to determine the range of evaluation of the density function. But
## this will fail for distributions with bigger absolute shape
## value. Instead I will use the range of the supplied data
plot.range <- seq( x.lim.rev[ 1 ] - bin.width,
x.lim.rev[ 2 ] + x$par[ 1 ]* 10, 0.01 )
plot.data <- data.frame(
x.plot = -1* plot.range,
y.plot = gpd.density( x$par, x$threshold * -1,
plot.range ) )
} else if ( x$control$extreme.type == "max" ){
## In the beginning I used the coefficients of the GP distribution
## to determine the range of evaluation of the density function. But
## this will fail for distributions with bigger absolute shape
## value. Instead I will use the range of the supplied data
plot.range <- seq( x.lim[ 1 ] - bin.width,
x.lim[ 2 ] + x$par[ 1 ]* 10, 0.01 )
plot.data <- data.frame(
x.plot = plot.range,
y.plot = gpd.density( x$par, x$threshold,
plot.range ) )
} else {
stop( "Unexpected extreme.type in provided object" )
}
## Removing all NaN to have less points and warnings
plot.data <- plot.data[
!is.nan( plot.data$y.plot ), ]
## Cutting of the density at values below the threshold in the tails
## Lower tail
plot.data <- plot.data[
( plot.data$x.plot >
plot.data$x.plot[ which.max( plot.data$y.plot ) ]
) | ( plot.data$y.plot > threshold.pdf.plot ), ]
## Upper tail
plot.data <- plot.data[
( plot.data$x.plot <
plot.data$x.plot[ which.max( plot.data$y.plot ) ]
) | ( plot.data$y.plot > threshold.pdf.plot ), ]
## If the first and last entry is not a zero, add one!
## Else the density won't look like a density at all.
if ( plot.data$y.plot[ 1 ] > 0 ){
plot.data <- data.frame(
x.plot = c( 2* plot.data$x.plot[ 1 ] -
plot.data$x.plot[ 2 ],
plot.data$x.plot ),
y.plot = c( 0, plot.data$y.plot ) )
}
if ( plot.data$y.plot[ nrow( plot.data ) ] > 0 ){
plot.data[ nrow( plot.data ) + 1, ] <-
c( 2* plot.data$x.plot[ nrow( plot.data ) ] -
plot.data$x.plot[ nrow( plot.data ) - 1 ],
0 )
}
plot.lim <- c( min( plot.data$x.plot ),
max( plot.data$x.plot ) )
## Using the threshold value defined above the range of
## the density plot will be determined
plot.lim <- c( plot.data[[ 1 ]][
which.min( abs( plot.data[[ 2 ]][
1 : which.max( plot.data[[ 2 ]] ) ] -
threshold.pdf.plot ) ) ],
plot.data[[ 1 ]][ which.min( (
plot.data[[ 2 ]][ which.max( plot.data[[ 2 ]] )
: length( plot.data[[ 2 ]] ) ] -
threshold.pdf.plot ) ) +
which.max( plot.data[[ 2 ]] ) - 1 ] )
## Determining the limits of the actual plot
if ( plot.lim[ 1 ] > x.lim[ 1 ] )
plot.lim[ 1 ] <- x.lim[ 1 ] - abs( x.lim[ 1 ] )* 0.05
if ( plot.lim[ 2 ] < x.lim[ 2 ] )
plot.lim[ 2 ] <- x.lim[ 2 ] + abs( x.lim[ 2 ] )* 0.05
plot.lim <- c( min( plot.lim[ 1 ],
min( x.data ) - bin.width ),
max( plot.lim[ 2 ],
max( x.data) + bin.width ) )
## Whenever the density gets bigger than one, cut it
if ( max( plot.data$y.plot ) > 1 ){
y.lim <- c( 0, ( max( stats::density( x.data )$y )* 1.5 ) )
} else {
y.lim <- NULL
}
ggplot() + geom_histogram(
data = data.frame( x = x.data ),
colour = grDevices::rgb( .098, .098, .44 ),
alpha = 1,
bins = bin.factor* length( x.data ),
aes( x = x.data, y = ..density..,
fill = "#7171EC" ) ) +
geom_polygon( data = plot.data, alpha = 0.7,
colour = grDevices::rgb( .098, .098, .44 ),
aes( x = x.plot, y = y.plot,
fill = grDevices::rgb( 1, .55, 0 ) ) ) +
scale_fill_manual( values = c( "#7171EC",
grDevices::rgb( 1, .55, 0 ) ),
labels = c( "Histogram", "Fitted GP" ) ) +
theme_bw() + ylab( "Density" ) +
coord_cartesian( xlim = plot.lim, ylim = y.lim ) +
theme( legend.title = element_blank() )
return( last_plot() )
}
## End of plot.R
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.