Implement mortalityTable.observed, add Austrian population mortality as example

DESCRIPTION

@@ -7,6 +7,7 @@ Authors@R: c(person("Reinhold", "Kainhofer", role=c("aut", "cre"), email="reinho
 Author: Reinhold Kainhofer [aut, cre]
 Maintainer: Reinhold Kainhofer <reinhold@kainhofer.com>
 URL: https://gitlab.open-tools.net/R/r-mortality-tables
+Encoding: UTF-8
 Depends:
     ggplot2,
     methods,
@@ -24,24 +25,27 @@ Description: Classes to implement and plot cohort life tables
     merged life tables.
 License: GPL (>= 2)
 RoxygenNote: 6.1.0.9000
-Collate: 'mortalityTable.R'
+Collate: 
+    'mortalityTable.R'
     'mortalityTable.period.R'
     'mortalityTable.ageShift.R'
     'ageShift.R'
-    'mortalityTable.observed.R'
-    'mortalityTable.joined.R'
+    'fillAges.R'
     'mortalityTable.mixed.R'
+    'mortalityTable.improvementFactors.R'
+    'mortalityTable.trendProjection.R'
+    'deathProbabilities.R'
+    'periodDeathProbabilities.R'
+    'mortalityTable.joined.R'
+    'getOmega.R'
+    'pensionTable.R'
+    'utilityFunctions.R'
+    'mortalityTable.observed.R'
     'ages.R'
     'baseTable.R'
     'baseYear.R'
-    'fillAges.R'
-    'pensionTable.R'
     'commutationNumbers.R'
-    'mortalityTable.improvementFactors.R'
-    'mortalityTable.trendProjection.R'
-    'deathProbabilities.R'
     'getCohortTable.R'
-    'getOmega.R'
     'getPeriodTable.R'
     'lifeTable.R'
     'makeQxDataFrame.R'
@@ -50,7 +54,6 @@ Collate: 'mortalityTable.R'
     'mortalityTable.MakehamGompertz.R'
     'mortalityTable.Weibull.R'
     'mortalityTable.deMoivre.R'
-    'periodDeathProbabilities.R'
     'mortalityTable.jointLives.R'
     'mortalityTables.list.R'
     'mortalityTables.load.R'
@@ -61,6 +64,5 @@ Collate: 'mortalityTable.R'
     'setLoading.R'
     'setModification.R'
     'undampenTrend.R'
-    'utilityFunctions.R'
     'whittaker.mortalityTable.R'
 VignetteBuilder: knitr

NAMESPACE

@@ -26,6 +26,7 @@ export(mortalityTable.deMoivre)
 export(mortalityTable.improvementFactors)
 export(mortalityTable.jointLives)
 export(mortalityTable.mixed)
+export(mortalityTable.observed)
 export(mortalityTable.once)
 export(mortalityTable.onceAndFuture)
 export(mortalityTable.period)
@@ -33,7 +34,9 @@ export(mortalityTable.trendProjection)
 export(mortalityTable.zeroes)
 export(mortalityTables.list)
 export(mortalityTables.load)
+export(pT.calculateTotalMortality)
 export(pT.getSubTable)
+export(pT.recalculateTotalMortality)
 export(pT.setDimInfo)
 export(pensionTable)
 export(pensionTables.list)
@@ -51,6 +54,7 @@ exportClasses(mortalityTable.deMoivre)
 exportClasses(mortalityTable.improvementFactors)
 exportClasses(mortalityTable.jointLives)
 exportClasses(mortalityTable.mixed)
+exportClasses(mortalityTable.observed)
 exportClasses(mortalityTable.period)
 exportClasses(mortalityTable.trendProjection)
 exportClasses(pensionTable)

R/fillAges.R

@@ -17,6 +17,7 @@ fillAges = function(probs = c(), givenAges = c(), neededAges = NULL, fill = NA_r
         result = rep(fill, length(neededAges))
         providedAges = intersect(neededAges, givenAges)
         result[match(providedAges, neededAges)] = probs[match(providedAges, givenAges)]
+        setNames(result, neededAges)
         result
     } else {
         probs

R/getOmega.R

@@ -30,9 +30,3 @@ setMethod("getOmega", "mortalityTable.mixed",
 #           function(object) {
 #               getOmega(object@table1);
 #           })
-
-# #' @describeIn getOmega Return the maximum age of the joined life table
-# setMethod("getOmega", "mortalityTable.observed",
-#           function(object) {
-#               max(object@ages, na.rm = TRUE);
-#           })

R/mortalityTable.observed.R

-#' @include mortalityTable.R
+#' @include mortalityTable.R utilityFunctions.R getOmega.R periodDeathProbabilities.R deathProbabilities.R ages.R
 NULL
 
-# #' Class mortalityTable.observed - Life table from actual observations
-# #'
-# #' A cohort life table described by actual observations (data frame of PODs
-# #' per year and age)
-# #'
-# #' @slot data    The observations
-# #' @slot years   The observation years
-# #' @slot ages    The observation ages
-# #'
-# #' @export mortalityTable.observed
-# #' @exportClass mortalityTable.observed
-# mortalityTable.observed = setClass(
-#     "mortalityTable.observed",
-#     slots = list(
-#         data = "data.frame",
-#         years = "numeric",
-#         ages = "numeric"
-#     ),
-#     prototype = list(
-#         data = data.frame(),
-#         years = c(),
-#         ages = c()
-#     ),
-#     contains = "mortalityTable"
-# )
-# asdf
+#' Class mortalityTable.observed - Life table from actual observations
+#'
+#' A cohort life table described by actual observations (data frame of PODs
+#' per year and age)
+#'
+#' @slot data    The observations
+#' @slot years   The observation years
+#' @slot ages    The observation ages
+#'
+#' @export mortalityTable.observed
+#' @exportClass mortalityTable.observed
+mortalityTable.observed = setClass(
+ "mortalityTable.observed",
+ slots = list(
+     deathProbs = "data.frame",
+     years = "numeric",
+     ages = "numeric"
+ ),
+ prototype = list(
+     deathProbs = data.frame(),
+     years = c(),
+     ages = c()
+ ),
+ contains = "mortalityTable"
+)
+
+#' @describeIn ages Return the ages of the life table
+setMethod("ages", "mortalityTable.observed",
+          function(object) {
+              object@ages
+          })
+
+#' @describeIn getOmega Return the maximum age of the life table
+setMethod("getOmega", "mortalityTable.observed",
+          function(object) {
+              max(object@ages, na.rm = TRUE);
+          })
+
+#' @describeIn mT.round Return the life table with the values rounded to the given number of digits
+setMethod("mT.round", "mortalityTable.observed",
+          function(object, digits = 8) {
+              o = callNextMethod()
+              o@data = round(o@data, digits = digits)
+              o
+          })
+
+
+
+# Solution to convert vector of integer values to string, with ranges for subsequent numbers:
+# https://stackoverflow.com/questions/16911773/collapse-runs-of-consecutive-numbers-to-ranges
+findIntRuns <- function(run) {
+    rundiff <- c(1, diff(run))
+    difflist <- split(run, cumsum(rundiff != 1))
+    runs = unlist(lapply(difflist, function(x) {
+        if (length(x) %in% 1:2) as.character(x) else paste0(x[1], "-", x[length(x)])
+    }), use.names = FALSE)
+    paste0(runs, collapse = ",")
+}
+
+#' @describeIn periodDeathProbabilities Return the (period) death probabilities
+#'             of the life table for a given observation year
+#'             If the observed mortality table does not provide data
+#'             for the desired period, the period closest to the
+#'             `Period` argument will be used and a warning printed.
+setMethod("periodDeathProbabilities", "mortalityTable.observed",
+          function(object, ..., ages = NULL, Period = 1975) {
+              if (is.null(ages)) {
+                  ages = ages(object)
+              }
+              col = which.min(abs(object@years - Period))
+              if (object@years[col] != Period) {
+                  warning("periodDeathProbabilities: Desired Period ", Period,
+                          " of observed mortalityTable not available, using closest period ",
+                          object@years[[col]], ".\nAvailable periods: ", findIntRuns(object@years))
+              }
+
+
+              # find the given year that is closest to the desired year:
+              #
+              fillAges(
+                  object@modification(object@deathProbs[,col] * (1 + object@loading)),
+                  givenAges = ages(object),
+                  neededAges = ages)
+          })
+
+#' @describeIn deathProbabilities Return the (cohort) death probabilities of the
+#'                                life table given the birth year (if needed)
+setMethod("deathProbabilities","mortalityTable.observed",
+          function(object,  ..., ages = NULL, YOB = 1975) {
+              if (is.null(ages)) {
+                  ages = ages(object);
+              }
+              years = YOB + ages;
+              yearcols = sapply(years, function(y) which.min(abs(object@years - y)))
+              agerows = match(ages, object@ages)
+
+              ## Check if all desired years are available
+              if (sum(abs(object@years[yearcols] - years)) > 0) {
+                  warning("deathProbabilities: Not all observation years ", findIntRuns(years),
+                          " of observed mortalityTable are available, using closest observations.\nAvailable periods: ", findIntRuns(object@years))
+              }
+
+              qx = object@deathProbs[cbind(agerows, yearcols)] * (1 + object@loading);
+              fillAges(object@modification(qx), givenAges = ages(object), neededAges = ages)
+          })
+
+
+#'@describeIn mT.cleanup Clean up the internal data of the mortality table
+setMethod("mT.cleanup", "mortalityTable.observed",
+          function(object) {
+              o = callNextMethod()
+              o@ages = unname(o@ages)
+              o@deathProbs = unname(o@deathProbs)
+              o@years = unname(o@years)
+              o
+          })

data-raw/Austria_Population_Observed_prepare.R

+library(tidyverse)
+library(openxlsx)
+filename = file.path("data-raw", "Austria_Population_Observed_StatistikAustria.xlsx")
+wb = openxlsx::loadWorkbook(filename)
+
+loadSheet = function(wb, sheet = "2017") {
+    if (as.numeric(sheet) >= 2002) {
+        startRow = 8
+        cols = c(1,2,8,14)
+        colNames = c("Alter", "M", "F", "U")
+    } else {
+        startRow = 13
+        cols = c(1,2,8)
+        colNames = c("Alter", "M", "F")
+    }
+
+    data = readWorkbook(wb, sheet = sheet, startRow = startRow, colNames = FALSE, rowNames = FALSE, cols = cols) %>%
+        `colnames<-`(colNames) %>%
+        filter(!is.na(M), !is.na(F)) %>%
+        mutate(Alter = as.integer(Alter), Jahr = as.integer(sheet)) %>%
+        gather(key = Geschlecht, value = qx, -Alter, -Jahr) %>%
+        select(Jahr, Alter, Geschlecht, qx) %>%
+        as.tibble
+
+    data
+}
+
+
+qx.AT_Pop_observed = bind_rows(sapply(sheets(wb), loadSheet, wb = wb, simplify = FALSE))
+
+for (g in c("M", "F", "U")) {
+    qx.AT_Pop_observed %>%
+        filter(Geschlecht == g) %>%
+        acast(Alter ~ Jahr, value.var = "qx") %>%
+        write.csv(file = file.path("inst", "extdata", paste0("Austria_Population_Observation_", g, ".csv")))
+}
+
+

inst/extdata/MortalityTables_Austria_PopulationObserved.R

+stopifnot(require(methods), require(utils), require(MortalityTables)) # MortalityTable classes; new; Excel reader
+
+
+###############################################################################
+### Gesamtbevölkerung Österreich: Bevölkerungsprognose bis 2080 (mittleres Szenario)
+### Datenquelle: Statistik Austria
+###############################################################################
+
+
+AT.pop.obs.M = utils::read.csv(system.file("extdata", "Austria_Population_Observation_M.csv", package = "MortalityTables"), check.names = FALSE, row.names = 1);
+AT.pop.obs.F = utils::read.csv(system.file("extdata", "Austria_Population_Observation_F.csv", package = "MortalityTables"), check.names = FALSE, row.names = 1);
+AT.pop.obs.U = utils::read.csv(system.file("extdata", "Austria_Population_Observation_U.csv", package = "MortalityTables"), check.names = FALSE, row.names = 1);
+
+mort.AT.observed.male = mortalityTable.observed(
+    name = "Österreich Männer Beobachtung",
+    deathProbs = AT.pop.obs.M,
+    ages = as.integer(rownames(AT.pop.obs.M)),
+    years = as.integer(colnames(AT.pop.obs.M)),
+    data = list(
+        dim = list(sex = "m", collar = "Gesamtbevölkerung", type = "Beobachtung", data = "official", year = "1947-2017")
+    )
+)
+mort.AT.observed.female = mortalityTable.observed(
+    name = "Österreich Frauen Beobachtung",
+    deathProbs = AT.pop.obs.F,
+    ages = as.integer(rownames(AT.pop.obs.F)),
+    years = as.integer(colnames(AT.pop.obs.F)),
+    data = list(
+        dim = list(sex = "f", collar = "Gesamtbevölkerung", type = "Beobachtung", data = "official", year = "1947-2017")
+    )
+)
+mort.AT.observed.unisex = mortalityTable.observed(
+    name = "Österreich Unisex Beobachtung",
+    deathProbs = AT.pop.obs.U,
+    ages = as.integer(rownames(AT.pop.obs.U)),
+    years = as.integer(colnames(AT.pop.obs.U)),
+    data = list(
+        dim = list(sex = "u", collar = "Gesamtbevölkerung", type = "Beobachtung", data = "official", year = "1947-2017")
+    )
+)
+
+rm(AT.pop.obs.M, AT.pop.obs.F, AT.pop.obs.U)
+
+###############################################################################
+
+# mortalityTables.load("Austria*")
+# plot(mort.AT.forecast.male, mort.AT.forecast.female, AVOe1996R.male, AVOe2005R.male, AVOe1996R.female, AVOe2005R.female, YOB = 2000)
+# plotMortalityTrend(mort.AT.forecast.male, mort.AT.forecast.female, AVOe1996R.male, AVOe2005R.male, AVOe1996R.female, AVOe2005R.female, Period = 2002)

man/mortalityTable.observed-class.Rd

+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/mortalityTable.observed.R
+\docType{class}
+\name{mortalityTable.observed-class}
+\alias{mortalityTable.observed-class}
+\alias{mortalityTable.observed}
+\title{Class mortalityTable.observed - Life table from actual observations}
+\description{
+A cohort life table described by actual observations (data frame of PODs
+per year and age)
+}
+\section{Slots}{
+
+\describe{
+\item{\code{data}}{The observations}
+
+\item{\code{years}}{The observation years}
+
+\item{\code{ages}}{The observation ages}
+}}
+