using-the-mortalityTables-package.Rmd 18.6 KB
Newer Older
1
---
2
title: "Using the MortalityTables Package"
3
4
5
author: "Reinhold Kainhofer, reinhold@kainhofer.com"
date: "`r Sys.Date()`"
output: 
6
    rmarkdown::html_vignette: 
7
8
9
10
11
        toc: true
        toc_depth: 3
        fig_width: 7
        fig_height: 5
vignette: >
12
  %\VignetteIndexEntry{MortalityTables}
13
14
15
16
17
18
19
20
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r echo = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
```

21
The MortalityTables package provides the `mortalityTable` base class and
22
some derived classes to handle diffferent types of mortality tables (also 
23
called life tables), mainly
24
25
26
27
28
29
30
used for life insurance. Additionally it provides a plot function to compare
multiple life tables either directly using the absolute mortalities in
log-linear plots or using relative mortalities as percentages of a given
reference table.

## Types of Life Tables

31
Provided types of mortality tables are:
32
33

* Base class
34
    : Class `mortalityTable`
35
* Period life table
36
    : Class `mortalityTable.period (ages, deathProbs, ..., baseYear=2000)`
37
38
39
    : Death probabilities observed / predicted for one observation year;
      No dependency on the bith year is assumed.
* Cohort life table using age-specific trends
40
    : Class `mortalityTable.trendProjection`
41
42
43
    : Death probabilities of a given base year are projected into the future
      using age-specific trends $\lambda_x$. The death probability of an $x$-year old in year
      `baseYear + n` is calculated as:
Reinhold Kainhofer's avatar
Reinhold Kainhofer committed
44
          $$q_x^{(baseYear+n)} = q_x^{(baseYear)} \cdot e^{-n\cdot\lambda_x}$$
45
    : Consequently, the death probabilities for a person born in year `YOB` can be calculated as
Reinhold Kainhofer's avatar
Reinhold Kainhofer committed
46
        $$q_x^{YOB} = q_x^{(base)} \cdot e^{-(YOB+x-baseYear)\cdot \lambda_x}$$
47
* Cohort life table approximation using age shift
48
    : Class `mortalityTable.ageShift`
49
50
    : Death probabilities for cohort $YOB$ are obtained by using death probabilities
      for cohort $X$ and modifying the technical age with a birth-year dependent shift:
Reinhold Kainhofer's avatar
Reinhold Kainhofer committed
51
52
          $$q_x^{YOB} = q_{x+shift(YOB)}^{(base)}$$
<!-- * Observed life table -->
53
<!--     : Class `mortalityTable.observed` -->
Reinhold Kainhofer's avatar
Reinhold Kainhofer committed
54
55
<!--     : Death probabilities observed during several years. The probabilities are -->
<!--       stored as a matrix with observation year and age as dimensions. -->
56
* Mixed life table
57
    : Class `mortalityTable.mixed`
58
59
60
61
    : Arithmetic mean of two life tables with given weights. This approach is
      often used to generate unisex life tables by mixing male and female
      mortalities with given weights (e.g. 70:30 or 40:60)
* Cohort life table using age-specific improvement factors
62
    : Class `mortalityTable.improvementFactors`
63
    : Project base life table using age-specific improvement factors.
64
65
66
67
* Pension tables
    : Class `pensionTable`
    : Transition probabilities for a four-state pension model (active, invalid, 
      retirement and death, with a possible widow in the event of death).
68

69
## Loading the MortalityTables package
70
```{r message=FALSE}
71
library("MortalityTables")
72
73
74
75
76
77
```

## Provided Data Sets

The package provides several real-life life tables published by census bureaus 
and actuarial associations around the world. You can use the function 
78
`mortalityTables.list` to list all available datasets (if no argument is given)
79
or all datasets that match the given pattern (wildcard character is *). You can 
80
then use `mortalityTables.load` to load either one single data set or all 
Reinhold Kainhofer's avatar
Reinhold Kainhofer committed
81
datasets that match the pattern.
82
83
84

```{r}
# list all available data sets
85
mortalityTables.list()
86
87

# list all datasets for Austria
88
mortalityTables.list("Austria_*")
89
90

# Load the German annuity table DAV 2004-R
91
mortalityTables.load("Germany_Annuities_DAV2004R")
92
93

# Load all Austrian data sets
94
mortalityTables.load("Austria_*")
95
96
97
98
99
100
101
102
103
104
105
```


In the next few sections we will always use some of the provided life tables 
for demonstration purposes. 


## Working with life table objects

### Plotting life tables

106
The package provides several functions to plot lifetables:
107

108
* `plotMortalityTables(table1, table2, ...)`
109
    : A log-linear plot comparing all given life tables.
110
* `plotMortalityTableComparisons(table1, table2, ..., reference=reftable)`
111
    : Plot the given life tables as percentages relative to the reference table
112
113
114
* `plotMortalityTrend(table1, table2, ..., YOB, Period)`
    : Plot the yearly mortality improvement factors (for either the given 
      observation year `Period` or the birth-year `YOB`)
115
    
116
These functionalities are also combined into the S3 plot function for the 
117
mortalityTable class, so you can usually just call plot on the mortality tables.
118
119
If the `trend = TRUE` argument is given, `plotMortalityTrend` is used,
if the `reference` argument is given, `plotMortalityTableComparisons` is used, 
120
otherwise `plotMortalityTables` is called.
121
122
123
124
```{r}
# Log-linear plot comparing some Austrian census tables
plot(mort.AT.census.1951.male, mort.AT.census.1991.male, 
     mort.AT.census.2001.male, mort.AT.census.2011.male, 
125
     legend.position = c(1,0))
126
127
128
129

# Relative death probabilities in percentage of the latest census
plot(mort.AT.census.1951.male, mort.AT.census.1991.male, 
     mort.AT.census.2001.male, 
130
     reference = mort.AT.census.2011.male, legend.position = c(1,0.75), ylim = c(0,4))
131
132
133
134
135
136
137
138
```

For cohort life tables, the plot functions also take either the `YOB` or the 
`Period` parameter to plot either the cohort death probabilities for the given 
birth year or the period death probabilities for the given observation year.

```{r}
# Comparison of two Austrian annuity tables for birth year 1977
139
plot(AVOe1996R.male, AVOe2005R.male, YOB = 1977, title = "Comparison for YOB=1977")
140
141

# Comparison of two Austrian annuity tables for observation year 2020
142
plot(AVOe1996R.male, AVOe2005R.male, Period = 2020, title = "Comparison for observation year 2020")
143
144
145
146
147
148
149
150
151
152
153
154

```

### Obtaining period and cohort death probabilities

To obtain death probabilities from all the different types of tables, there are two functions:
    
* `deathProbabilities`: Returns the (cohort) death probabilities of the life table given the birth year
* `periodDeathProbabilities`: Returns the (period) death probabilities of the life table for a given
    observation year

```{r message=FALSE}
155
mortalityTables.load("Austria_Annuities")
156
# Get the cohort death probabilities for Austrian Annuitants born in 1977:
157
qx.coh1977 = deathProbabilities(AVOe2005R.male, YOB = 1977)
158
159

# Get the period death probabilities for Austrian Annuitants observed in the year 2020:
160
qx.per2020 = periodDeathProbabilities(AVOe2005R.male, Period = 2020)
161
162
163
164
165
166
```

These functions return the death probabilities as a simple, numeric R vector. 

There are two similar functions that return the death probabilities as a period life table object that can be used with all other functions provided by this package:

167
168
* `getCohortTable`: Get a `mortalityTable` object describing the death probabilities for people born in the given year
* `getPeriodTable`: Get a `mortalityTable` object describing the death probabilities observed in the given year
169
170

```{r}
171
# Get the cohort death probabilities for Austrian Annuitants born in 1977 as a mortalityTable.period object:
172
table.coh1977 = getCohortTable(AVOe2005R.male, YOB = 1977)
173
174

# Get the period death probabilities for Austrian Annuitants observed in the year 2020:
175
table.per2020 = getPeriodTable(AVOe2005R.male, Period = 2020)
176
177

# Compare those two in a plot:
178
plot(table.coh1977, table.per2020, title = "Comparison of cohort 1977 with Period 2020", legend.position = c(1,0))
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217

```

Not surprisingly, at 43 years the two death probabilities cross, because in 2020
the person born 1977 is 43 years old, so the $q_x$ refer to the same person. 
However, for younger ages, the period 2020 probabilities are lower, because 
the mortality improvement for those younger ages has much less time in the 
cohort 1977 table. For ages above 43 the cohort table describes the mortality 
further into the future than 2020, so there is more improvement and thus lower 
death probabilities for the cohort life table.




### Other data extraction functions from life tables

| function               | description |
|:---------------------- |:---------------------------------------------------|
|`ages(table)`           | Returns the vector of ages, for which the life table can provide death probabilities |
|`getOmega(table)`       | Returns the maximum age, for which the life table can provide dath probabilities |
|`ageShift(table, YOB)`  | Returns the age shift for the given year of birth |
|`baseTable(table)`      | Returns the base table, from which the table projects (for cohort tables) |
|`baseYear(table)`       | Returns the year of the base table |
|`lifetable(table, YOB)`  | Returns the cohort death probabilities as a `lifetable` object for use with the lifecontingencies package|





## Creating a life table object

### Period life tables
Period death probabilities are the simplest type of life table, giving the 
probabilities of death observed during the
corresponding year (the "period"). The death probabilities of different ages
refer to different persons, being of the corresponding ages in the observation
year. All that is needed to create a period life table are the death probabilities
and the corresponding ages:
```{r}
218
219
lt = mortalityTable.period(name = "Sample period lifetable", ages = 1:99, deathProbs = exp(-(99:1)/10))
plot(lt, title = "Simple log-linear period mortality table")
220
221
222
223
224
deathProbabilities(lt)

```


Reinhold Kainhofer's avatar
Reinhold Kainhofer committed
225
226
<!-- ### Observed life tables -->
<!-- The observations for the given years -->
227

Reinhold Kainhofer's avatar
Reinhold Kainhofer committed
228
<!-- TODO -->
229
230
231
232


### Cohort life tables with trend projection

Reinhold Kainhofer's avatar
Reinhold Kainhofer committed
233
234
235
236
237
238
239
240
A cohort life table with trend projection needs the following parameters:

* The base table $q_x^{(base)}$ (death probabilities) for the given base period as a vector
* Age-specific trend factors $\lambda_x$ as a vector
* The base year (numeric)
* 

```{r}
241
atPlus2 = mortalityTable.trendProjection(
Reinhold Kainhofer's avatar
Reinhold Kainhofer committed
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
    name = "Austrian Census Males 2011, 2% yearly trend",
    baseYear = 2011,
    deathProbs = deathProbabilities(mort.AT.census.2011.male),
    ages = ages(mort.AT.census.2011.male),
    trend = rep(0.02, length(ages(mort.AT.census.2011.male)))
)
```

Some life tables do not assume a constant age-specific trend over time, but rather
assume that the currently observed high mortality improvements are just a 
temporary effect, so the current trend is in effect only for some time and 
then reduces to some kind of long-term trend.

There are two conceptual approaches: One is to use a trend dampening function
that is simply applied to the starting trend. So, while the initial trend might 
be 3\%, i.e. the projection will use `(ObservationYear-BaseYear) * OriginalYear`, 
over time it will assume the value 
`dampeningFunction(ObservationYear-BaseYear) * OriginalTrend`. The dampening 
function in this case gives the cumulated trend effect from the base year until 
the observation year.
262
To implement this trend reduction with the MortalityTables package, simply pass
Reinhold Kainhofer's avatar
Reinhold Kainhofer committed
263
264
265
a one-argument function as the `dampingFunction` slot to the class, the argument 
will be the number of years from the base year (NOT the calendar year!):
```{r}
266
atPlus2.damp = mortalityTable.trendProjection(
Reinhold Kainhofer's avatar
Reinhold Kainhofer committed
267
268
269
270
271
272
273
274
    name = "Austrian M '11, 2% yearly, damping until 2111",
    baseYear = 2011,
    deathProbs = deathProbabilities(mort.AT.census.2011.male),
    ages = ages(mort.AT.census.2011.male),
    trend = rep(0.02, length(ages(mort.AT.census.2011.male))),
    # damping function: 2011: full effect, linear reduction until yearly trend=0 in 2111:
    # 2011: 100%, 2012: 99%, 2013: 98% => For 2013 we have a cumulative trend 
    # of 297% instead of 300% for three full yearly trends!
275
    dampingFunction = function(n) { n - n * (n + 1) / 2 / 100 }
Reinhold Kainhofer's avatar
Reinhold Kainhofer committed
276
277
)

278
plot(mort.AT.census.2011.male, atPlus2, atPlus2.damp, YOB = 2011, legend.position = c(0.8,0.75))
Reinhold Kainhofer's avatar
Reinhold Kainhofer committed
279
280
281
282
283
284
285
286
287
288
289
290
```

The other approach is to assume that instead of the initial trend, after some 
time a second trend (slot trend2) takes over. In this case, the `dampingFunction`
slot is again a one-argument function that now gives the weight of the first trend, while `1-dampingFunction(year)` will give the weight of the second trend. As the weights 
will be applied for the whole period from the base- to the observation year, the weights
need to be cumulated and normalized. 

The argument
in this case is the actual calendar year (not the year since the base year like it was in the one-trend case above!)

```{r}
291
atPlus2.damp2 = mortalityTable.trendProjection(
Reinhold Kainhofer's avatar
Reinhold Kainhofer committed
292
293
294
295
296
297
298
299
300
301
    name = "Austrian M '11, 2% yearly, 1% long-term",
    baseYear = 2011,
    deathProbs = deathProbabilities(mort.AT.census.2011.male),
    ages = ages(mort.AT.census.2011.male),
    trend = rep(0.02, length(ages(mort.AT.census.2011.male))),
    trend2 = rep(0.01, length(ages(mort.AT.census.2011.male))),
    # damping function interpolates between the two trends: 
    # until 2021 trend 1, from 2031 trend 2, linearly beteen
    dampingFunction = function(year) { 
        if (year <= 2021) 1
302
303
        else if (year > 2031) 14.5/(year - 2011)
        else 1 - (year - 2021)*(year - 2021 + 1) / 20 / (year - 2011)
Reinhold Kainhofer's avatar
Reinhold Kainhofer committed
304
305
306
    }
)

307
plot(mort.AT.census.2011.male, atPlus2, atPlus2.damp, atPlus2.damp2, YOB = 2011, legend.position = c(0.8,0.75))
Reinhold Kainhofer's avatar
Reinhold Kainhofer committed
308
```
309
310
311

### Cohort life tables with age-shift

Reinhold Kainhofer's avatar
Reinhold Kainhofer committed
312
313
314
315
316
317
318
319
320
321
322
323
324
325
Age-shifted cohort life tables are an approximation to full cohort life tables.
Full cohort life tables apply a trend or improvment factors to the death 
probabilities of a base year to obtail death probabilities for a given birth year.
Age-shifting rather modifies the age of the corresponding person and uses the 
same, unmodified base table for all cohorts. Basically, it works like this:

> A 60-year old born in 1950 has the same death probability as a 50-year old 
> born in 1900, so instead of looking at the cohort 1950, we can look at the 
> cohort 1900 and for a person born 1950 we treat him as if he were 10 years 
> younger.

So, an age-shifted cohort life table just needs the base table and for each 
birth year the amount the age is modified.

326
327
328
329
330
331
332
333
334
335
For those people, who think visually, age shifting works on the death 
probabilities as following: A normal trend moves the $q_x$ curve downwards. 
Age-shifting approximates this by shifting the $q_x$ curve to the right without
modifying its values.

The following example clearly shows this, with the blue curve being the base 
table for YOB 2011. A full trend projection moves the curve down to the green line,
while age-shifting moves the base curve to the right so that it coincides as 
much as possible with the exact (green) line.

Reinhold Kainhofer's avatar
Reinhold Kainhofer committed
336
```{r}
337
baseTableShift = getCohortTable(atPlus2, YOB = 2011);
Reinhold Kainhofer's avatar
Reinhold Kainhofer committed
338
339
baseTableShift@name = "Base table of the shift (YOB 2011)"

340
atShifted = mortalityTable.ageShift(
Reinhold Kainhofer's avatar
Reinhold Kainhofer committed
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
    name = "Approximation with age shift",
    baseYear = 2011,
    deathProbs = deathProbabilities(baseTableShift),
    ages = ages(baseTableShift),
    ageShifts = data.frame(
        shifts = c(
            rep( 0, 3), 
            rep(-1, 3), 
            rep(-2, 3), 
            rep(-3, 3), 
            rep(-4, 3), 
            rep(-5, 3), 
            rep(-6, 3)
        ),
        row.names = 2011:2031
    )
)

359
ageShift(atShifted, YOB = 2021)
Reinhold Kainhofer's avatar
Reinhold Kainhofer committed
360

361
plot(baseTableShift, atPlus2, atShifted, YOB = 2021, legend.position = c(0.8,0.75))
Reinhold Kainhofer's avatar
Reinhold Kainhofer committed
362
363
364
365
366
367
368
369
```

As one can see, for ages above 40 years, the table with 2% yearly trend and the
corresponding age-shifted table have roughly the same mortalities. Below 40 years, 
the two are very different, so this approximation through age-shifting should 
really be used with extreme care!


370
371
372
373
374
375
376
377
378


## Modifying life table objects

### Copying life tables

Life tables are simple pass-by-value S4 objects, so copying works by simple assignment. 

```{r}
379
b = AVOe2005R.female 
380
381
382
b@name = "Modified Copy"
# only b is modified, not the original table
b@modification = function(qx) pmax(qx, 0.01)  
383
plot(AVOe2005R.female, b, YOB = 2000)
384
385
```

386
### Adding a security loading to the raw probabilities
387

388
389
390
391
392
When calculating premiums for life insurance contracts, one often needs to add 
a certain security loading on the raw death probabilities (e.g. 10% increased
death probabilities) to account for statistical fluctuations. This can be easily
done with the `setLoading` function that returns a copy of the given table and 
adds the given security loading.
393
394

```{r}
395
396
397
AVOe2005R.female.sec = setLoading(AVOe2005R.female, loading = 0.1);
# Make sure the modified table has a new name, otherwise plots might break
AVOe2005R.female.sec@name = "Table with 10% loading"
398
plot(AVOe2005R.female, AVOe2005R.female.sec, title = "Original and modified table")
399
400
401
402
403
```

### Adding a modification to the raw probabilities

Some uses require post-processing of the death probabilities, like adding a lower
404
bound for the death probabilities. To achive this, all `mortalityTable`-derived 
405
406
407
408
classes have a slot `modification` that takes a function that is passed the vector 
of death probabilities.

```{r}
409
AVOe2005R.female.mod = setModification(AVOe2005R.female, modification = function(qx) pmax(0.03, qx));
410
411
# Make sure the modified table has a new name, otherwise plots might break
AVOe2005R.female.mod@name = "Modified table (lower bound of 3%)"
412
plot(AVOe2005R.female, AVOe2005R.female.mod, title = "Original and modified table")
413
414
```

415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
## Pension Tables
Pension tables generalize mortality tables in that the state space is increased 
from two states (alive / dead) to four states (active / invalidity or realy 
retirement / old age retirement / dead). As a consequence, there is no longer 
just one transition probability, but multiple.

Possible states are:
* active: healty, no pension, typically paying some kin of premium
* incapacity: disablity pension (in most cases permanent), not working, early pension
* retirement: old age pension, usually starting with a fixed age
* dead
  * Widow/widower pension (if a widow exists at the time of death)

Correspondingly, the `pensionTable` class offers the following slots describing 
transition probabilities for the corresponding state transitions (by a 
`mortalityTable`-derived object):
* `qxaa`:  death probability of actives (active -> dead)
* `ix`:    invalidity probability (active -> incapacity)
* `qix`:   death probability of invalid (invalid -> dead)
* `rx`:    reactivation probability (incapacity -> active)
* `apx`:   retirement probability (active -> retirement), typically 1 for a fixed age
* `apx`:   retirement probability of invalids (invalid -> retirement), typically 0 or 1 for a fixed age
* `qpx`:   death probability of retired (retired -> dead)
* `hx`:    probability of a widow at moment of death (dead -> widow), y(x) age differene
* `qxw`:   death probability of widows/widowers
* `qgx`:   death probability of total group (irrespective of state)

All functions that handle `mortalityTable` object can be used with these slots.

Additionally, the following functions are provided to obtain the set of all 
transition probabilities in one data frame:
* `transitionProbabilities(pension_table, YOB)`
* `periodTransitionProbabilities(pension_table, Period)`