R/plotMortalityTrend.R

In MortalityTables: A Framework for Various Types of Mortality / Life Tables

#' Plot the trends of multiple mortality tables (life tables) in one chart
#'
#' \code{plotMortalityTrend} prints the trends of multiple life tables (objects
#' of child classes of \code{mortalityTable}) in one plot, with a legend showing
#' the names of the tables.
#'
#' @param data First life table to be plotted. Either a \code{data.frame} generated by \code{makeQxDataFrame} or a \code{mortalityTable} object
#' @param ... Additional life tables to be plotted (if \code{data} is a \code{mortalityTable} object)
#' @param aes Optional aesthetics to append or override the default. The default aesthetics will always be applied first and provide defaults for x, y and color. This argument can be used to override the defaults or append other aesthetics.
#' @param xlim X-axis limitatation (as a two-element vector)
#' @param ylim Y-axis limitatation (as a two-element vector)
#' @param xlab X-axis label (default: "Alter")
#' @param ylab Y-axis label (default: "Sterbewahrscheinlichkeit q_x relativ zu ....")
#' @param title The plot title
#' @param legend.position The position of the legend (default is \code{c(0.9,0.1)})
#' @param legend.justification The justification of the legend (default is \code{c(1,)})
#' @param legend.key.width The keywith of the lines in the  legend (default is \code{unit(25,"mm")})
#' @param legend.title Title of the legend (\code{NULL} to hide)
#' @param ages Plot only the given ages
#'
#' @examples
#' # Load the Austrian aunnity data
#' mortalityTables.load("Austria_Annuities")
#'
#' # Compare the trends of these tables
#' plotMortalityTrend(AVOe2005R.male, AVOe2005R.female, AVOe1996R.male, AVOe1996R.female,
#'     Period = 2002, title = "Trends of Austrian Annuity Tables")
#' # For tables with a non-constant trend, the Period and YOB can be used to compare
#' # the age-specific trends that apply to the death probabilities during a given
#' # period or for a given birth year
#' plotMortalityTrend(AVOe2005R.male, AVOe2005R.female, AVOe1996R.male, AVOe1996R.female,
#'     YOB = 1950, title = "Trends of Austrian Annuity Tables for cohort YOB=1950")
#' plotMortalityTrend(AVOe2005R.male, AVOe2005R.female, AVOe1996R.male, AVOe1996R.female,
#'     YOB = 2000, title = "Trends of Austrian Annuity Tables for cohort YOB=2000")
#' plotMortalityTrend(AVOe2005R.male, AVOe2005R.female, AVOe1996R.male, AVOe1996R.female,
#'    Period = 1999, title = "Trends of Austrian Annuity Tables for Period 2002")
#' plotMortalityTrend(AVOe2005R.male, AVOe2005R.female, AVOe1996R.male, AVOe1996R.female,
#'    Period = 2030, title = "Trends of Austrian Annuity Tables for Period 2030")
#'#' @import scales
#' @export
plotMortalityTrend = function(
    data, ...,
    aes = NULL,
    ages = NULL,
    xlim=NULL, ylim=NULL,
    xlab=NULL, ylab=NULL,
    title = "",
    legend.position=c(0.9, 0.9), legend.justification = c(1, 1), legend.title = "Sterbetafel",
    legend.key.width = unit(25, "mm")
) {

    if (!is.data.frame(data)) {
        data = makeMortalityTrendDataFrame(data, ...);
    }
    if (!is.null(ages)) {
        data = data[data$x %in% ages,]
    }
    if (missing(xlab)) xlab = "Alter";
    if (missing(ylab)) ylab = expression(paste("Sterblichkeitstrend ", lambda[x]));

    pl = ggplot(data, aes(x = x, y = y, color = group))
    if (!is.null(aes)) {
        pl = pl + aes
    }
    pl = pl +
        theme_bw() +
        theme(
            plot.title = element_text(size = 18, face = "bold"),
            legend.title = element_text(size = 14, face = "bold.italic"),
            # legend in bottom right corner of the plot
            legend.justification = legend.justification, legend.position = legend.position,
            # No box around legend entries
            legend.key = element_blank(),
            legend.key.width = legend.key.width,
            legend.background = element_rect(colour = "gray50", linetype = "solid")
        ) +
        geom_line() +
        scale_y_continuous(
            name = ylab,
            #      breaks = scales::trans_breaks('log10', function(x) 10^x),
            #      labels = scales::trans_format('log10', scales::math_format(10^.x))
            #minor_breaks = log(c(sapply(x, function(x) seq(0, x, x/10))), 10)
        ) +
        scale_x_continuous(
            name = xlab,
            #breaks = function (limits) scales::trans_breaks('', function(x) 10^x),
            breaks = function(limits) seq(max(min(limits), 0), max(limits), 5),
            minor_breaks = function(limits) seq(max(round(min(limits)), 0), round(max(limits)), 1)
            #labels = scales::trans_format('log10', scales::math_format(10^.x))
        ) +
        coord_cartesian(xlim = xlim, ylim = ylim) +
        xlab("Alter") + labs(colour = legend.title);
    if (title != "") {
        pl = pl + ggtitle(title);
    }
    pl
}

makeMortalityTrendDataFrame = function(..., YOB = 1972, Period = NULL) {
    # If reference is given, normalize all probabilities by that table!
    data = unlist(list(...));
    if (is.null(data)) return(data.frame(x = double(), y = double(), group = character()))
    # names(data) = lapply(data, function(t) t@name);

    if (missing(Period) || is.null(Period)) {
        data = lapply(data, function(t) {
            if (is.data.frame(t@data$dim) || is.list(t@data$dim)) {
                data.frame(x = ages(t), y = `names<-`(mortalityImprovement(t, YOB = YOB), NULL), group = t@name, as.data.frame(t@data$dim))
            } else {
                data.frame(x = ages(t), y = `names<-`(mortalityImprovement(t, YOB = YOB), NULL), group = t@name)
            }
        });
    } else {
        data = lapply(data, function(t) {
            if (is.data.frame(t@data$dim) || is.list(t@data$dim)) {
                data.frame(x = ages(t), y = `names<-`(mortalityImprovement(t, Period = Period), NULL), group = t@name, as.data.frame(t@data$dim))
            } else {
                data.frame(x = ages(t), y = `names<-`(mortalityImprovement(t, Period = Period), NULL), group = t@name)
            }
        });
    }

    names(data) = NULL
    data <- as.data.frame(do.call("rbind.expand", data))
    data
}


globalVariables(c("x", "y", ".x"))


# mortalityTables.load("Austria_*")
# plotMortalityTrend(AVOe1996R.male, AVOe1996R.female, AVOe2005R.male, AVOe2005R.female, YOB=1972, title="Austrian Annuity Tables, YOB=1972 (for cohort tables)")
#
# plotMortalityTables(mort.AT.census.2001.male, AVOe2005R.male, AVOe2005R.female, YOB=1972, title="Comparison Austrian Tables")
#  plotMortalityTables(getCohortTable(AVOe2005R.male, YOB=1972), getCohortTable(AVOe2005R.male, YOB=2016), title="Comparison Austrian Tables")

# plotMortalityTrend(EttlPagler.male@qx, AVOe1999P.male@qx, AVOe2008P.male@qx, YOB = 2003, title="Sterblichkeitstrends der Pagler-Tafeln")