From 3222a0d0a9c5a1d70b9eec342e5ba735113a523f Mon Sep 17 00:00:00 2001
From: Reinhold Kainhofer <reinhold@kainhofer.com>
Date: Thu, 25 Aug 2016 00:35:00 +0200
Subject: [PATCH] Add plotValuationTableComparisons function
Basically the same as plotValuationTables, except that all mortality rates are normalized by a given mortality table and the y-axis is linear (not logarithmic). This allows easy comparison of tables.
---
R/plotValuationTables.R | 136 +++++++++++++++++++++++++++++++++++-----
1 file changed, 120 insertions(+), 16 deletions(-)
diff --git a/R/plotValuationTables.R b/R/plotValuationTables.R
index e056851..1e05048 100644
--- a/R/plotValuationTables.R
+++ b/R/plotValuationTables.R
@@ -1,20 +1,67 @@
#' @export
-makeQxDataFrame = function(..., YOB=1972, Period=NA) {
- data=list(...);
- names(data) = lapply(data, function(t) t@name);
- if (missing(Period)) {
- cat("Year of birth: ", YOB, "\n");
- data = lapply(data, function(t) cbind(x=ages(t), y=deathProbabilities(t, YOB=YOB)))
- } else {
- cat("Period: ", Period,"\n");
- data = lapply(data, function(t) cbind(x=ages(t), y=periodDeathProbabilities(t, Period=Period)))
- }
+#' If reference is given, normalize all probabilities by that table!
+makeQxDataFrame = function(..., YOB=1972, Period=NA, reference=NULL) {
+ data=list(...);
+ names(data) = lapply(data, function(t) t@name);
+ reference_ages = NULL;
+
+ if (missing(Period)) {
+ cat("Year of birth: ", YOB, "\n");
+ if (!missing(reference)) {
+ reference_ages = ages(reference);
+ reference = deathProbabilities(reference, YOB=YOB);
+ }
+ data = lapply(data, function(t) {
+ str(deathProbabilities(t, YOB=YOB));
+ normalize_deathProbabilities(
+ cbind(x=ages(t), y=deathProbabilities(t, YOB=YOB)),
+ reference = reference,
+ referenceAges = reference_ages)
+ });
+ } else {
+ cat("Period: ", Period,"\n");
+ if (!missing(reference)) {
+ reference_ages = ages(reference);
+ reference = periodDeathProbabilities(reference, Period=Period);
+ }
+ data = lapply(data, function(t) {
+ str(deathProbabilities(t, YOB=YOB));
+ normalize_deathProbabilities(
+ cbind(x=ages(t), y=periodDeathProbabilities(t, Period=Period)),
+ reference = reference,
+ referenceAges = reference_ages)
+ });
+ }
+
+ list.names = names(data)
+ lns <- sapply(data, nrow)
+ data <- as.data.frame(do.call("rbind", data))
+ data$group <- rep(list.names, lns)
+ data
+}
+
+normalize_deathProbabilities = function (data, reference=NULL, referenceAges=NULL) {
+ if (missing(reference) || missing(referenceAges) || is.null(reference) || is.null(referenceAges)) {
+ return(data);
+ }
+# browser();
+ # Find which ages exist in both and obtain those indices from the data and the reference list:
+ useages = intersect(data[,"x"], referenceAges)
+ dataindices = match(useages, data[,"x"])
+ refindices = match(useages, referenceAges)
+
+ # Find which ages in data do NOT exist in the reference ages (and are thus NOT normalized at all)
+ # Print a warning!
+ missingrefs = setdiff(data[,"x"], referenceAges)
+ if (length(missingrefs)>0) {
+ warning("Reference mortality table does not contain ages ",
+ missingrefs,
+ " required for normalization. These ages will not be normalized!")
+ }
- list.names = names(data)
- lns <- sapply(data, nrow)
- data <- as.data.frame(do.call("rbind", data))
- data$group <- rep(list.names, lns)
- data
+ # Now divide the data by the corresponding entries from the reference list
+ data[dataindices, "y"] = data[dataindices, "y"] / reference[refindices]
+ data
}
#' Plot multiple valuation tables (life tables) in one plot
@@ -28,7 +75,7 @@ makeQxDataFrame = function(..., YOB=1972, Period=NA) {
#' @param legend.key.width The keywith of the lines in the legend (default is \code{unit(25,"mm")})
#'
#' @export
-plotValuationTables = function(data, ..., title = "", legend.position=c(0.9,0.1), legend.key.width = unit(25, "mm")) {
+plotValuationTables = function(data, ..., xlim=NULL, ylim=NULL, title = "", legend.position=c(0.9,0.1), legend.key.width = unit(25, "mm")) {
if (!is.data.frame(data)) {
data = makeQxDataFrame(data, ...);
}
@@ -60,6 +107,7 @@ plotValuationTables = function(data, ..., title = "", legend.position=c(0.9,0.1)
#labels = scales::trans_format('log10', scales::math_format(10^.x))
) +
+ coord_cartesian(xlim=xlim, ylim=ylim) +
annotation_logticks(sides="lr") +
xlab("Alter") + labs(colour="Sterbetafel");
if (title != "") {
@@ -67,6 +115,62 @@ plotValuationTables = function(data, ..., title = "", legend.position=c(0.9,0.1)
}
pl
}
+
+plotValuationTableComparisons = function(data, ..., xlim=NULL, ylim=NULL, title = "", legend.position=c(0.9,0.1), legend.key.width = unit(25, "mm"), reference=NULL) {
+ # If no reference mortality table is given, use the first table (data if its a valuation table)
+ if (missing(reference)) {
+ if (inherits(data, "valuationTable")) {
+ reference = data;
+ } else {
+ reference = NULL;# TODO;
+ }
+ }
+ refname=reference@name;
+ yAxisLabel = expression(paste("Sterbewahrscheinlichkeit ", q[x], " relativ zu ", refname));
+yAxisLabel[1][3]="asfd"
+yAxisLabel
+ if (!is.data.frame(data)) {
+ data = makeQxDataFrame(data, ..., reference=reference);
+ }
+
+ pl = ggplot(data, aes(x = x, y = y, colour = data$group)) +
+ 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=c(1,0), 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() +
+ coord_cartesian(xlim=xlim, ylim=ylim) +
+ scale_y_continuous(
+# #name=substitute(paste("Sterbewahrscheinlichkeit ", q[x], " relativ zu ", refname), env=list(refname=refname))#,
+ name=substitute(refname^3, env=list(refname=refname)),
+ labels=percent
+# # 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="Alter",
+ #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))
+
+ ) +
+ # annotation_logticks(sides="lr") +
+ xlab("Alter") + labs(colour="Sterbetafel");
+ if (title != "") {
+ pl = pl + ggtitle(title);
+ }
+ pl
+}
+
#
# plotValuationTables(mort.AT.census.1869.male, mort.AT.census.1869.female, mort.AT.census.2011.male, mort.AT.census.2011.female, AVOe2005R.male, AVOe2005R.female, YOB=1972,title="Vergleich österreichische Sterbetafeln, YOB=1972 (bei Generationentafeln)")
#
--
GitLab