What is lemur?

lemur (“Life expectancy monitor upscaled in R”) is an R package and Shiny application for scenario analysis of mortality. Given the cause-of-death distribution of a population, you ask “what if?” questions – what would life expectancy be if cardiovascular mortality fell by 50%? – and the package recomputes the life table, shows you the gains and losses at every age, and decomposes the change by age and cause of death. It also compares cause-of-death profiles and life tables between regions, sexes and time periods.

Everything runs on the bundled Global Burden of Disease 2023 (GBD 2023) estimates from IHME: 216 regions from Inf to -Inf (the calendar-year-2023 data included), three sex categories (male, female and the two sexes combined) and 18 broad cause groups, arranged in abridged life tables of 22 age groups.

The tool is hosted by the HASS Digital Research Hub (HDRH) at the Australian National University (ANU).

Installation

Once you have an R session open (in RStudio or the R console), install the package straight from GitHub:

lemur requires R version 4.3.0 or newer.

# 1. Install the pak package, if you do not have it yet
install.packages("pak")

# 2. Install lemur from GitHub (latest version on the main branch)
pak::pak("mpascariu/lemur")

The package bundles the GBD 2023 datasets (exposed as data_gbd_lt(), data_gbd_cod() and data_gbd_sdg()), covering deaths and life tables from 1990 to 2023, so you can run the examples on the help pages and the analysis functions without a database connection.

The data

The three accessor functions return pre-factorized data.tables:

• data_gbd_lt() – one abridged life table per region, sex and year, with the classic columns mx, qx, ax, lx, dx, Lx, Tx and ex.
• data_gbd_cod() – death counts by region, sex, year, age group and cause (deaths column).
• data_gbd_sdg() – the same, grouped into the cause categories used for tracking the UN Sustainable Development Goals.

Let’s look at Romania in 2021 (both sexes combined). An excerpt of the life table – survivors lx, deaths dx, probability of dying qx and remaining life expectancy ex – at selected ages:

library(lemur)

# helper used throughout this document, defined at first use
slice_table <- function(tbl, region, sex = "both", year = 2021) {
  filter(tbl, region == .env$region, sex == .env$sex, period == .env$year)
}

lt_ro  <- slice_table(data_gbd_lt(), "Romania")
cod_ro <- slice_table(data_gbd_cod(), "Romania")
cardio <- c("Ischemic Heart Disease", "Stroke", "Other Cardiovascular")
total_deaths <- sum(cod_ro$deaths)

Life expectancy at birth in Romania in 2021 was 72.81 years, with 75,235 of every 100,000 newborns surviving to age 65. The corresponding cause-of-death table counts 334,246.4 deaths in that year, distributed over 18 causes and 22 age groups.

The largest single block – COVID-19 – accounts for 21.7% of all deaths. (2021 was a pandemic year, so COVID-19 tops the list.) This is the kind of target a scenario analysis will aim at.

Workflow I – “What if” scenario analysis

The core idea: mortality at age x from cause c can be reduced or increased by a factor, and the all-cause life table is rebuilt from the modified cause-specific risks. This is what the Shiny app does every time you move a slider, and the same functions are exported for use in scripts.

The pipeline has three steps:

1. build_cod_matrix() turns the long cause-of-death table into an age × cause matrix of death shares (each row sums to 1).
2. modify_cod_table() (or modify_life_table()) applies a percentage change – as a scalar, a vector, or a full age × cause matrix.
3. The modified risks are converted back into a life table.

Let’s model a realistic question: what if cardiovascular mortality in Romania in 2021 were halved, at every age? In the GBD classification, cardiovascular deaths are split across three cause groups – ischaemic heart disease, stroke and other cardiovascular – together worth 43.7% of all deaths. The change is a matrix of zeros with -50 in the three cardiovascular columns:

causes <- as.character(unique(cod_ro$cause_name))
ages   <- unique(lt_ro$x)

M <- matrix(0, nrow = length(ages), ncol = length(causes),
            dimnames = list(ages, causes))
M[, cardio] <- -50

scenario <- prepare_data(
  cod = cod_ro, lt = lt_ro,
  region1 = "Romania", region2 = "Romania", sex = "both",
  cod_change = M, mode = "cod"
)
lt_after  <- scenario$lt_final
cod_after <- scenario$cod_final

prepare_data() is the exact function the app runs: it filters the inputs, applies the change and returns the initial and final life tables and cause-of-death tables. The result:

Life expectancy at birth rises by 2.04 years (72.8 → 74.8), and roughly 73,000 fewer deaths occur in a year. Note that the deaths drop by about half of the cardiovascular deaths – the other causes are untouched.

Where does this gain in life expectancy come from? plot_change() plots the remaining life expectancy of the first life table minus the second at every age: a marker to the right of the dashed zero line is a gain, a marker to the left a loss.

plot_change(L1 = lt_after, L2 = lt_ro,
            title = "Halving cardiovascular mortality: gains by age",
            subtitle = "Romania, 2021, both sexes")

The gains are strikingly uniform: halving cardiovascular mortality adds 2.04 years of remaining life expectancy at birth, 2.05 at age 30, 1.87 at 60 and still 1.4 years at 80, before tapering to at age 100. Because the reduction applies proportionally at every age, the saved person-years accrue to the whole population rather than piling up in old age; the gain fades only where little life remains to be gained. The cause-of-death profile after the change tells the other half of the story:

plot_cod(cod_after, perc = TRUE)

Compare this with the initial distribution: the cardiovascular share of all deaths falls from 43.7% to 27.9%. Every other cause keeps its death count exactly – only the three cardiovascular columns changed. That is the “what if” contract: a small set of levers moved, everything else held constant.

Finally, we can decompose the 2.04-year gain into age-by-cause contributions with decompose_by_cod() and draw the result with plot_decompose():

dec_scenario <- decompose_by_cod(lt_ro, lt_after, cod_ro, cod_after)
plot_decompose(dec_scenario, by = "both")

As expected, the entire gain is attributed to the three modified causes. Positive values favor the second population passed – here the post-change table – which is why the whole figure lies above the zero line. The stack shows the age pattern: the age group contributing the most is [75,80) (0.29 years), the ages where cardiovascular deaths are concentrated. Each bar is the sum of the per-cause contributions, and the sum of the whole figure equals the total gain to the last decimal.

Workflow II – decomposing differences between populations

The same toolkit answers comparative questions: why does male life expectancy differ between Romania and Mexico in 2021? First, the gap:

lt1  <- slice_table(data_gbd_lt(), "Romania", sex = "male")
lt2  <- slice_table(data_gbd_lt(), "Mexico",  sex = "male")
cod1 <- slice_table(data_gbd_cod(), "Romania", sex = "male")
cod2 <- slice_table(data_gbd_cod(), "Mexico",  sex = "male")

Male life expectancy at birth is 69.24 years in Romania and 67.38 years in Mexico – a gap of 1.86 years, favoring Romania. decompose_by_age() splits this gap into age contributions (using the direct-indirect method of Andreev, Shkolnikov and Begun, 2002):

dec_age <- decompose_by_age(lt1, lt2)

The sign convention is simple: positive values favor the second population passed to the function – here Mexico – and negative ones favor the first (Romania). The age pattern shows where the advantage is won and where it is lost: the [60,65) age group contributes the most in Mexico’s favor (0.272 years), while the [0,1) group works the other way (-0.544 years, favoring Romania). The column sums exactly to the total gap.

decompose_by_cod() refines this by cause, using stepwise replacement: it replaces the mortality of one population with the other’s, cause by cause and age by age, and reads off the life-expectancy effect of each replacement. The symmetrical = TRUE (default) and direction = "up" arguments follow the recommendations of Andreev et al. (2002), averaging the two replacement directions.

dec <- decompose_by_cod(lt1, lt2, cod1, cod2)

In absolute terms the largest single contribution is COVID-19 (-1.67 years, in Romania’s favor): Mexico’s male COVID-19 mortality in 2021 far exceeds Romania’s. The largest factor pulling the other way is Stroke (0.87 years, in Mexico’s favor) – Romanian males’ stroke mortality sits well above Mexico’s. The full picture – each cause at each age – is best seen in the stacked plot:

plot_decompose(dec, by = "both")

Reading the bars: segments above the zero line favor the second population (Mexico), segments below favor the first (Romania). The colours follow the package-wide cause palette, so the same cause carries the same colour in every figure. Summing every segment recovers the total gap to the last decimal – the decomposition is exact, not an approximation.

Cause-of-death differences can also be viewed directly with plot_cod() in comparison mode, or collapsed to a single horizontal bar with plot_decompose(by = "cod"):

plot_decompose(dec, by = "cod")

The Shiny application

Everything above is wrapped in an interactive dashboard, launched with a single call:

lemur::run_app()          # local data mode (bundled GBD data)

The interactive dashboard, launched with lemur::run_app().

The dashboard offers five analysis modes:

In every mode, moving a slider or ticking a cause re-runs the pipeline shown above; a 250 ms debounce batches input changes into a single consistent update, cached reactives avoid recomputing unchanged queries, and all charts render natively in plotly (roughly 0.1 s per interaction). A “Data” tab exposes the underlying tables – initial and final life tables, death counts, the decomposition and the reduction matrix – and the “Methods” tab explains the algorithms. The map tab shows the selected region on an interactive leaflet map with 2021 population, life expectancy and fertility indicators.

The SDG modes use the bundled data_gbd_sdg() table, which groups the same death counts into the SDG cause categories:

This makes it easy to ask, for instance, how much life expectancy would improve if the SDG 3.4 target – a one-third reduction of premature mortality from non-communicable diseases – were achieved in a given region.

References

• Every exported function has a worked example: ?modify_life_table, ?decompose_by_cod, ?plot_decompose, ?data_gbd_lt.
• The decomposition methods follow Andreev, Shkolnikov and Begun (2002), Algorithm for decomposition of differences between aggregate demographic measures, Demographic Research 7, 499-522.
• The data are the Global Burden of Disease Study 2023 (GBD 2023) results, Institute for Health Metrics and Evaluation (IHME) (https://vizhub.healthdata.org/gbd-results/). The calendar-year-2023 data is included; deaths and life tables span 1990-2023.
• Source code, issues and releases: https://github.com/mpascariu/lemur.