%\VignetteIndexEntry{inedemogR: A Comprehensive Tutorial} %\VignetteEngine{R.rsp::tex} %\VignetteEncoding{UTF-8} \documentclass[11pt]{article} \usepackage[T1]{fontenc} \usepackage[utf8]{inputenc} \usepackage{mathptmx} % Times-compatible text and math font \usepackage[margin=1in]{geometry} \usepackage{amsmath} \usepackage{amssymb} \usepackage{booktabs} \usepackage{longtable} \usepackage{array} \usepackage{xcolor} \usepackage{listings} \usepackage{enumitem} \usepackage{caption} \usepackage[colorlinks=true,linkcolor=blue,citecolor=blue,urlcolor=blue]{hyperref} \usepackage{parskip} \definecolor{codebg}{rgb}{0.96,0.96,0.96} \definecolor{codekw}{rgb}{0.0,0.35,0.55} \definecolor{codestr}{rgb}{0.55,0.0,0.0} \definecolor{codecom}{rgb}{0.4,0.4,0.4} \lstdefinelanguage{Rlang}{ morekeywords={function,if,else,for,while,repeat,break,next,return,TRUE,FALSE,NULL,NA,library}, sensitive=true, morecomment=[l]{\#}, morestring=[b]", morestring=[b]' } \lstset{ language=Rlang, basicstyle=\ttfamily\footnotesize, keywordstyle=\color{codekw}\bfseries, commentstyle=\color{codecom}\itshape, stringstyle=\color{codestr}, backgroundcolor=\color{codebg}, frame=single, framerule=0.3pt, rulecolor=\color{gray!40}, breaklines=true, breakatwhitespace=false, columns=flexible, showstringspaces=false, xleftmargin=2pt, xrightmargin=2pt, aboveskip=8pt, belowskip=8pt } \newcommand{\pkg}[1]{\textbf{#1}} \newcommand{\fn}[1]{\texttt{#1}} \title{\pkg{inedemogR}: Tidy Access to Spanish INE Demographic Data\\[4pt] \large A Comprehensive Tutorial} \author{J.\ R.\ Caro-Barrera} \date{\today} \begin{document} \maketitle \tableofcontents \newpage \section{Introduction} \pkg{inedemogR} provides tidy, harmonized, reproducible access to demographic data published by the Spanish National Statistics Institute (Instituto Nacional de Estad\'{i}stica, INE), retrieved live via the official \pkg{ineapir} API wrapper, with spatial integration via \pkg{mapSpain} and plotting built on \pkg{ggplot2}. The package is organized around two complementary data systems that serve different purposes. They are deliberately kept parallel rather than merged into one interface, because they operate at genuinely different levels of demographic detail. \begin{description}[leftmargin=2.2cm,style=nextline] \item[System A] \fn{get\_ine\_demog()}, \fn{get\_ine\_geo()}, \fn{list\_ine\_indicators()}, \fn{plot\_ine\_map()}, \fn{map\_indicator()}, \fn{update\_ine\_data()}. Quick, multi-geography (municipality or province) totals for population, births, and deaths, with no age or sex breakdown. Best for national/regional choropleths and quick comparisons. \item[System B] \fn{get\_ine\_population()}, \fn{get\_ine\_births()}, \fn{get\_ine\_births\_by\_age()}, \fn{get\_ine\_deaths()}, \fn{compute\_exposure()}, \fn{compute\_death\_rates()}, \fn{build\_life\_tables()}, \fn{build\_abridged\_life\_tables()}, \fn{download\_ine\_data()}, and the summary indicators, decomposition, and charts built on top of them. Rigorous, province-level, age/sex-disaggregated analysis following the Human Mortality Database (HMD) Methods Protocol V6, as implemented for the Spanish subnational Human Mortality Database (SHMD) pipeline. This is the only system capable of producing dependency ratios, population pyramids, Lexis diagrams, period life tables, fertility schedules (ASFR/TFR/MAC/GRR/NRR), and mortality decomposition, since these all require an age breakdown that System A's data does not carry. \end{description} This tutorial is split into two parts, matching the two purposes a user of \pkg{inedemogR} typically has: \begin{itemize} \item \textbf{Part I (Section~\ref{sec:part1})} covers strictly data \emph{retrieval, cleaning, validation, and storage} --- the functions in \fn{get\_ine\_demog.R}, \fn{get\_ine\_geo.R}, \fn{population.R}, \fn{births.R}, \fn{deaths.R}, \fn{exposure.R}, \fn{death\_rates.R}, \fn{life\_tables.R}, \fn{abridged\_life\_tables.R}, \fn{download\_ine\_data.R}, and the caching layer in \fn{helpers.R}. \item \textbf{Part II (Section~\ref{sec:part2})} covers \emph{demographic analysis and visualization} --- the indicator functions in \fn{indicators.R}, the mortality-decomposition functions in \fn{decomposition.R}, and the charting/mapping functions in \fn{plot\_demog.R}, \fn{plot\_ine\_map.R}, and \fn{lexis.R} --- together with full worked examples. \end{itemize} Every mortality-pipeline calculation (exposure-to-risk, central death rates, life tables) is presented with the exact mathematical formulation implemented in the code, following the HMD Methods Protocol V6. \subsection{Installation} \begin{lstlisting} # From a local clone or source tarball: # install.packages("devtools") devtools::install() library(inedemogR) \end{lstlisting} \subsection{A note on live data} Every retrieval function in this package calls the real INE API (or, for \fn{get\_ine\_geo()}, the real \pkg{mapSpain} boundary service) over the network. There is no bundled offline dataset. All code listings in this tutorial are real, runnable calls against live INE data as of the 2023--2024 reference years used in the examples; results will naturally update as INE publishes new years. Section~\ref{sec:caching} explains how to avoid repeatedly re-downloading the same data during iterative analysis. \newpage \section{Part I --- Data Retrieval, Cleaning, and Storage} \label{sec:part1} \subsection{System A: quick multi-geography totals} \subsubsection{\fn{list\_ine\_indicators()}} Returns the internal registry (\fn{ine\_variables}) mapping each indicator code to the real INE table it is retrieved from, the geographic dimension column INE's API returns for that table, and the geographic granularity the table is actually published at. \begin{lstlisting} list_ine_indicators() #> # A tibble: 3 x 6 #> indicator name description id_table geo_var geo_level #> population_total Resident pop. ... 29005 Municipios municipality #> births_total Live births ... 6506 Provincias province #> deaths_total Deaths ... 6545 Provincias province \end{lstlisting} Migration indicators are not included in this release: INE's migration tables are origin$\times$destination$\times$year flow matrices, a different shape from the simple geography$\times$year stock tables backing the three indicators above, so wiring one in needs its own design pass rather than a registry one-liner. \fn{get\_ine\_demog()} still raises an explicit \fn{"Unknown indicator(s)"} error for any indicator code not in \fn{list\_ine\_indicators()}, migration or otherwise. \subsubsection{\fn{get\_ine\_demog()}} Retrieves one or more indicators via \fn{ineapir::get\_data\_table()} and tidies the result into one row per geography $\times$ year, one column per indicator. \begin{lstlisting} get_ine_demog(indicator, geo_level = NULL, year = NULL, region = NULL, sex = "Total", geometry = FALSE) \end{lstlisting} \begin{description}[leftmargin=2.4cm,style=nextline] \item[indicator] Character vector of indicator codes (see \fn{list\_ine\_indicators()}). \item[geo\_level] The geographic level to fetch; inferred from \fn{indicator} if \fn{NULL}. All requested indicators must share the same level (population is municipality-level; births/deaths are province-level only) --- requesting a mix across levels is a hard error. \item[year] Numeric/character year; \fn{NULL} returns only the latest published period. \item[region] Optional regex vector matched against the geography name. \item[sex] One of \fn{"Total"}, \fn{"Hombres"}, \fn{"Mujeres"}. \item[geometry] If \fn{TRUE}, joins \fn{get\_ine\_geo()}'s geometry and returns an \fn{sf} object. \end{description} \begin{lstlisting} # Population totals for every municipality, 2023 pop_totals <- get_ine_demog(indicator = "population_total", year = 2023) # Births and deaths together, province level, one call vitals <- get_ine_demog( indicator = c("births_total", "deaths_total"), year = 2023 ) \end{lstlisting} \subsubsection{\fn{get\_ine\_geo()}} Retrieves province or municipality boundaries via \pkg{mapSpain} (\fn{esp\_get\_prov()} / \fn{esp\_get\_munic()}), reprojected to ETRS89 / UTM zone 30N (EPSG:25830), with a \fn{GEOID} column that matches \fn{get\_ine\_demog()}'s output exactly so the two can be joined directly. \begin{lstlisting} get_ine_geo(geo_level = c("municipality", "province"), region = NULL, moveCAN = TRUE, can_gap_km = 60) \end{lstlisting} When \fn{moveCAN = TRUE} (the default), the Canary Islands are translated so they sit \fn{can\_gap\_km} kilometres from the nearest mainland coastline point (near C\'{a}diz/Huelva) rather than left in their true, far-southwest position --- the standard cartographic convention for compact national maps of Spain. The shifted islands' bounding box is attached as the \fn{"can\_box"} attribute, used by \fn{plot\_ine\_map()} and \fn{map\_indicator()} to draw an inset separator frame. \begin{lstlisting} prov_geo <- get_ine_geo(geo_level = "province") mun_geo <- get_ine_geo(geo_level = "municipality", region = "Sevilla") \end{lstlisting} \subsubsection{\fn{update\_ine\_data()}} Checks every indicator wired in the registry for its latest published year and refreshes a local \fn{.rds} cache (under \fn{tools::R\_user\_dir("inedemogR", "cache")} by default) only if INE has published a newer year than what is cached, or if \fn{force = TRUE}. \begin{lstlisting} update_ine_data(force = FALSE, data_dir = tools::R_user_dir("inedemogR", "cache")) \end{lstlisting} \subsection{System B: age/sex-disaggregated province data} \subsubsection{Caching layer (\fn{helpers.R})} \label{sec:caching} \fn{get\_ine\_population()}, \fn{get\_ine\_births()}, and \fn{get\_ine\_deaths()} share a common disk-cache wrapper, \fn{with\_ine\_cache()}. Fetching province-level age-specific deaths and population from INE's Tempus3 API (sometimes via a full CSV bulk-export fallback when the JSON API rejects an oversized query) is the slowest part of the package, often taking minutes; results are cached as \fn{.rds} files keyed by the call's arguments so repeated calls during iterative analysis are near-instant. Every System B retrieval function exposes the same three arguments: \begin{description}[leftmargin=2.2cm,style=nextline] \item[use\_cache] (default \fn{TRUE}) Return a previously cached result for the same arguments instead of re-fetching. \item[force] (default \fn{FALSE}) Ignore any existing cache and re-fetch from INE, refreshing the cache afterwards. \item[cache\_dir] Directory for the cache (default: \fn{tools::R\_user\_dir("inedemogR", "cache")}, the same location \fn{update\_ine\_data()} uses). \end{description} \begin{lstlisting} # First call: fetches live and caches the result. pop <- get_ine_population(n_periods = 10) # Second call, same arguments: returns instantly from cache. pop <- get_ine_population(n_periods = 10) # Force a fresh fetch, e.g. after INE has published a new year. pop <- get_ine_population(n_periods = 10, force = TRUE) \end{lstlisting} \subsubsection{\fn{get\_ine\_population()}} Retrieves Spanish January-1st provincial population by single-year age (0--100+) and sex from INE's Padr\'{o}n Municipal Continuo, per SHMD Protocol Section 12, Step 3. \begin{lstlisting} get_ine_population(table_id = 56945, n_periods = 30, use_cache = TRUE, force = FALSE, cache_dir = tools::R_user_dir("inedemogR", "cache")) \end{lstlisting} Returns \fn{list(data, qc)}: \fn{data} has columns \fn{ine\_code}, \fn{nuts3\_code}, \fn{nuts2\_code}, \fn{province\_name}, \fn{year}, \fn{age}, \fn{female}, \fn{male}, \fn{total}; \fn{qc} is the output of \fn{validate\_population()}, which checks for negative/missing counts, incomplete age ranges, duplicate rows, year gaps, and implausible age-to-age jumps ($>$50\%, informational only). \begin{lstlisting} pop <- get_ine_population() head(pop$data) pop$qc$passed \end{lstlisting} \subsubsection{\fn{get\_ine\_births()}} Retrieves annual live births by province and sex from INE's MNPN operation. \begin{lstlisting} get_ine_births(table_id = 6506, n_periods = 100, use_cache = TRUE, force = FALSE, cache_dir = tools::R_user_dir("inedemogR", "cache")) \end{lstlisting} Returns \fn{list(data, qc)}; \fn{data} has \fn{ine\_code}, \fn{nuts3\_code}, \fn{nuts2\_code}, \fn{province\_name}, \fn{year}, \fn{female}, \fn{male}, \fn{total}. \fn{validate\_births()} checks the same invariants as \fn{validate\_population()} (no age dimension here, so no age-range check). \subsubsection{\fn{get\_ine\_births\_by\_age()}} \label{sec:births-by-age} \fn{get\_ine\_births()} has no age-of-mother breakdown, so it cannot feed a true total fertility rate. \fn{get\_ine\_births\_by\_age()} retrieves the same MNPN operation's age-of-mother table instead (table 6508: province, single year of age of the mother 15--49 plus the two open intervals \fn{under15}/\fn{50plus}, and sex of the newborn) --- the input Section~\ref{sec:fertility-schedule}'s age-specific fertility rate schedule needs. \begin{lstlisting} get_ine_births_by_age(table_id = 6508, n_periods = 100, use_cache = TRUE, force = FALSE, cache_dir = tools::R_user_dir("inedemogR", "cache")) \end{lstlisting} Returns \fn{list(data, qc)}; \fn{data} has \fn{ine\_code}, \fn{nuts3\_code}, \fn{nuts2\_code}, \fn{province\_name}, \fn{year}, \fn{age} (integer, \fn{NA} for the \fn{under15} group), \fn{age\_group} (\fn{"under15"}, \fn{"15"} --- \fn{"49"}, \fn{"50plus"}), \fn{female}, \fn{male}, \fn{total}. INE's table omits a (province, age, sex-of-newborn, year) series entirely once its value is zero across the whole fetched window --- the same behavior already seen in \fn{compute\_death\_rates()} (Section~\ref{sec:death-rates}) --- so \fn{clean\_births\_by\_age()} builds the complete grid explicitly and fills absent combinations with zero rather than leaving them missing, before \fn{validate\_births\_by\_age()} checks for negative/missing counts, duplicate rows, and complete single-year age coverage (15--49) per province-year. \begin{lstlisting} births_age <- get_ine_births_by_age() births_age$data \end{lstlisting} \subsubsection{\fn{get\_ine\_deaths()}} Retrieves annual deaths by province and sex (INE's MNPD operation, table 6545) \emph{and} an age-specific companion series (table 6547), since the mortality pipeline (\fn{compute\_exposure()}, \fn{compute\_death\_rates()}) needs deaths by single-year age, $D(x,t)$, not just the yearly total $D(t)$. \begin{lstlisting} get_ine_deaths(table_id = 6545, age_table_id = 6547, n_periods = 100, use_cache = TRUE, force = FALSE, cache_dir = tools::R_user_dir("inedemogR", "cache")) \end{lstlisting} Returns \fn{list(data\_provinces, data\_national, qc, qc\_age)}: \fn{data\_provinces} is the age-specific series (\fn{nuts3\_code}, \fn{province\_name}, \fn{year}, \fn{age}, \fn{female}, \fn{male}, \fn{total}) that feeds the rest of the mortality pipeline; \fn{data\_national} is Spain's age-less annual total, used internally to cross-check the sum of provincial totals against the official national figure (\fn{validate\_deaths()}'s \fn{national\_vs\_provincial\_mismatch} check) and also validates the sex ratio at death lies in the plausible range $[0.90, 1.50]$. \begin{lstlisting} deaths <- get_ine_deaths() deaths$data_provinces # age-specific, feeds compute_exposure()/compute_death_rates() deaths$data_national # Spain total, used for QC cross-checks \end{lstlisting} \subsection{The mortality pipeline} The three functions below implement, in order, Steps 4--6 of the SHMD Protocol (equivalently, HMD Methods Protocol V6): exposure-to-risk, central death rates, and period life tables. Each stage consumes the previous stage's output directly. \subsubsection{Exposure-to-risk: \fn{compute\_exposure()}} Exposure-to-risk $E(x,t)$ approximates the person-years lived at age $x$ during calendar year $t$, combining the population stock at the start and end of the year with a Lexis-triangle correction for within-year deaths: \begin{equation} E(x,t) \;=\; \tfrac{1}{2}\bigl[P(x,t) + P(x,t+1)\bigr] \;+\; \tfrac{1}{2}\bigl[D^{U}(x,t) - D^{L}(x,t)\bigr] \end{equation} where $P(x,t)$ is the January~1 population aged $x$ in year $t$, and $D^{U}(x,t)$, $D^{L}(x,t)$ are deaths in the upper and lower Lexis triangles of the $(x,t)$ parallelogram --- i.e.\ deaths to the birth cohort $t-x$ (upper triangle, born before the start of year $t$, dying at age $x$) versus deaths to the birth cohort $t-x-1$ (lower triangle, born during year $t$ itself among those who reach age $x$). INE's Tempus3 death tables do not carry a birth-cohort split, so absent a \fn{cohort} column, \fn{compute\_exposure()} falls back to the documented HMD Methods Protocol V6 Appendix A approximation of an even 50/50 split ($D^{U}=D^{L}=D(x,t)/2$, which cancels the correction term entirely, i.e.\ reduces to a plain average of the two population stocks). \begin{lstlisting} compute_exposure(population, deaths) \end{lstlisting} \begin{description}[leftmargin=2.2cm,style=nextline] \item[population] \fn{get\_ine\_population()\$data}. \item[deaths] \fn{get\_ine\_deaths()\$data\_provinces}. \end{description} Because death registration lags population estimates by roughly a year, the most recent year in \fn{population} commonly has no matching rows in \fn{deaths} yet. That year is still used as the $P(x,t+1)$ boundary needed to compute the \emph{previous} year's exposure, but is not itself returned as an exposure row --- treating an entirely-absent year as zero deaths would mechanically produce zero mortality and, downstream, a nonsensical inflated life expectancy. \begin{lstlisting} pop <- get_ine_population() deaths <- get_ine_deaths() exposure <- compute_exposure(pop$data, deaths$data_provinces) exposure$qc$passed \end{lstlisting} \fn{validate\_exposure()} flags cells with no usable source population (informational), any remaining negative/NA exposure, and implausible exposure exceeding twice the source population (a generous plausibility bound). \subsubsection{Central death rates: \fn{compute\_death\_rates()}} \label{sec:death-rates} Given deaths $D(x,t)$ and exposure $E(x,t)$, the central death rate at single-year age $x$ in year $t$ is \begin{equation} m(x,t) \;=\; \frac{D(x,t)}{E(x,t)} \qquad \text{(1x1: single-year age, one-year period)} \end{equation} INE's age-specific deaths table only publishes rows for age$\times$year cells with at least one death; \fn{compute\_death\_rates()} starts from the complete $0$--\fn{MAX\_AGE} exposure grid and left-joins deaths, filling absent rows with zero deaths rather than silently dropping them --- a missing INE row means \emph{zero} deaths at that age, not a data gap, and dropping it would break the age-contiguity that the life-table recursion (Section~\ref{sec:lifetables}) requires. A pooled 5-year-age-group ($n=5$) rate is also computed: \begin{equation} m(x,n,t) \;=\; \frac{\sum_{k=0}^{4} D(x+k,t)}{\sum_{k=0}^{4} E(x+k,t)} \qquad \text{(5x1)} \end{equation} and an age-standardized death rate using the Eurostat European Standard Population 2013 (ESP2013) single-year age weights $w(x)$, with $\sum_x w(x) = 100{,}000$: \begin{equation} \mathrm{ASDR}(t) \;=\; \sum_{x} m(x,t)\, w(x) \end{equation} expressed in deaths per 100,000 of the ESP2013 standard population, allowing direct comparison of mortality levels across provinces or years with different underlying age structures. \begin{lstlisting} compute_death_rates(deaths, exposure) \end{lstlisting} Returns \fn{list(mx\_1x1, mx\_5x1, asdr, qc)}. \fn{validate\_death\_rates()} flags negative/NA/infinite rates and any rate $\geq 1$ below age 90 (implausible for a working-age population, informational). \begin{lstlisting} rates <- compute_death_rates(deaths$data_provinces, exposure$data) rates$mx_1x1 # age x year x province mortality schedule rates$asdr # age-standardized rate, one row per province-year \end{lstlisting} \subsubsection{Period life tables: \fn{build\_life\_tables()}} \label{sec:lifetables} \fn{build\_life\_tables()} constructs standard 1x1 period life tables for every province and year in \fn{mx\_1x1}, following HMD Methods Protocol V6 exactly. The construction proceeds in two stages: old-age smoothing, then the life-table recursion itself. \paragraph{Kannisto smoothing of old-age mortality.} Single-year mortality rates at the oldest ages are noisy (small populations, few deaths), so ages 80--110 are replaced with a smooth Kannisto model fitted on ages 80--95 and extrapolated to age 110: \begin{equation} \operatorname{logit} m(x) = \log\!\frac{m(x)}{1-m(x)} = \alpha + \beta x, \qquad x \in [80, 95] \end{equation} fitted by ordinary least squares, then evaluated for every age up to the open-ended terminal age (110): \begin{equation} m^{\text{Kannisto}}(x) = \frac{e^{\alpha + \beta x}}{1 + e^{\alpha + \beta x}} \end{equation} \paragraph{Andreev--Kingkade $a_0$.} The average person-years lived by those who die before age 1, $a_0$, uses the Andreev--Kingkade (2015) piecewise regression on the age-0 central death rate $m_0$ (replacing the older Coale--Demeny approximation), separately by sex: \begin{equation} a_0 = \begin{cases} 0.330 & \text{if } m_0 \geq 0.107 \\ 0.053 + 2.800\,m_0 & \text{if female and } m_0 < 0.107 \\ 0.045 + 2.684\,m_0 & \text{if male and } m_0 < 0.107 \end{cases} \end{equation} For the combined-sex table, $a_0$ is the mean of the female and male values. The remaining $a_x$ values follow the standard HMD convention: $a_1 = 0.4$ (age 1), $a_x = 0.5$ for ages $2$ and above. \paragraph{Life-table recursion.} Given $m_x$ and $a_x$ for every age $x = 0, 1, \dots, \omega$ (open interval at $\omega = 110$): \begin{align} q_x &= \frac{m_x}{1 + (1-a_x)m_x}, \qquad q_\omega = 1 \\ l_0 &= 100{,}000 \\ d_x &= l_x\, q_x \\ l_{x+1} &= l_x - d_x \\ L_x &= \begin{cases} l_{x+1} + a_x d_x & x < \omega \\ l_\omega / m_\omega & x = \omega \end{cases}\\ T_x &= \sum_{y \geq x} L_y \\ e_x &= T_x / l_x \end{align} The open-interval formula $L_\omega = l_\omega/m_\omega$ assumes $\omega$ is genuinely old (a high $m_\omega$); this is the reason the terminal age must be a realistic $\omega=110$ after Kannisto extrapolation rather than an artificially low cut-off, which would otherwise produce a nonsensical remaining-life-expectancy figure in the thousands of years. \begin{lstlisting} build_life_table(mx_df, sex = c("female", "male")) # single year, single sex build_life_tables(mx_1x1) # every province and year \end{lstlisting} \fn{build\_life\_tables()} returns \fn{list(fltper, mltper, bltper, qc)} --- female, male, and both-sex period life tables, each a tibble with \fn{nuts3\_code}, \fn{province\_name}, \fn{year}, \fn{age}, \fn{mx}, \fn{qx}, \fn{ax}, \fn{lx}, \fn{dx}, \fn{Lx}, \fn{Tx}, \fn{ex}. \fn{validate\_life\_table()} checks that $l_x$ is non-increasing in age, $L_x>0$, $e_x>0$, and $q_x \in [0,1]$. \begin{lstlisting} lt <- build_life_tables(rates$mx_1x1) lt$fltper[lt$fltper$age == 0, ] # e0 by province and year, female \end{lstlisting} As a live-verified example, this pipeline gives a period life expectancy at birth for A Coru\~{n}a, 2024, of $e_0 = 86.97$ years (female) and $81.21$ years (male) --- plausible, correct figures for a Spanish province in that year. \subsubsection{Abridged life tables: \fn{build\_abridged\_life\_tables()}} \label{sec:abridged} \fn{build\_abridged\_life\_tables()} builds standard abridged (5-year age group) period life tables directly from \fn{compute\_death\_rates()}'s \fn{mx\_5x1}, complementing \fn{build\_life\_tables()}'s single-year (1x1) tables --- useful for comparing against other agencies' published abridged tables, or when only grouped-age data is available. Age groups are \fn{"00-04"}, \fn{"05-09"}, \dots, \fn{"95-99"}, \fn{"100+"} (open interval). The one exception to building strictly from \fn{mx\_5x1} is the youngest group: INE's own age grouping merges age 0 with ages 1--4 into a single \fn{"00-04"} group, which would otherwise need an external, less-precise $_5a_0$ regression (e.g.\ Coale--Demeny) to handle the sharp within-interval infant-mortality gradient. Since \fn{mx\_1x1} already carries the exact single-year $m_0$ with its Andreev--Kingkade $a_0$ (Section~\ref{sec:lifetables}), \fn{build\_abridged\_life\_table()} instead builds an exact single-year sub-table for ages 0--4 and aggregates it, avoiding that approximation for the one interval where it matters most. For every other closed interval $[x, x+n)$ ($n=5$), the standard discrete abridged life-table identities apply (Preston, Heuveline \& Guillot, 2001, Ch.\ 3): \begin{align} {}_na_x &= n/2 \\ {}_nq_x &= 1 - e^{-n\, {}_nm_x} \end{align} The exponential form for ${}_nq_x$ (rather than the linear Chiang formula $n\, {}_nm_x / (1 + (n - {}_na_x)\, {}_nm_x)$ that \fn{build\_life\_table()} uses for single-year, $n=1$ intervals) is deliberate: the linear formula can exceed 1 once $n$ is as wide as 5 and mortality is high (old-age groups with ${}_nm_x \gtrsim 0.3$--$0.4$), which the exponential form never does. The remaining recursion matches Section~\ref{sec:lifetables} exactly, with the person-years-lived identity generalized to the interval width: \begin{equation} {}_nL_x = n\, l_{x+n} + {}_na_x\, {}_nd_x \qquad (x < \omega) \end{equation} \begin{lstlisting} build_abridged_life_table(mx_1x1, mx_5x1, sex = c("female", "male")) # one province/year build_abridged_life_tables(mx_1x1, mx_5x1) # every province/year \end{lstlisting} Returns \fn{list(fltper, mltper, bltper, qc)}, matching \fn{build\_life\_tables()}'s structure with \fn{age\_group}/\fn{age\_start}/ \fn{n} in place of \fn{age}. \fn{validate\_abridged\_life\_table()} checks the same invariants as \fn{validate\_life\_table()}. \begin{lstlisting} alt <- build_abridged_life_tables(rates$mx_1x1, rates$mx_5x1) alt$fltper[alt$fltper$age_group == "00-04", ] \end{lstlisting} As a live-verified cross-check, A Coru\~{n}a 2024's abridged $e_0$ (female) comes to $87.03$ years against the exact single-year figure of $86.97$ --- a $0.06$-year difference, consistent with the coarser 5-year grouping. \subsection{Exporting to files: \fn{download\_ine\_data()}} For users who prefer to work outside R (a spreadsheet, another statistical package, or HMD-compatible mortality software), \fn{download\_ine\_data()} runs any subset of the six pipeline stages and writes the results to a folder as CSV and/or HMD-format \fn{.txt} files. \begin{lstlisting} download_ine_data(out_dir, stages = c("births", "deaths", "population", "exposure", "mx", "life_tables"), format = c("csv", "txt"), n_periods = 30) \end{lstlisting} Later stages automatically pull in their dependencies (\fn{exposure} needs \fn{population} and \fn{deaths}; \fn{mx} needs \fn{deaths} and \fn{exposure}; \fn{life\_tables} needs \fn{mx}), since they are fetched along the way regardless. \begin{lstlisting} # Everything, both file formats, into ./ine_data result <- download_ine_data("ine_data") result$files # every path written result$life_tables$fltper # also available in memory, no re-read needed # Just births and deaths, CSV only download_ine_data("ine_data", stages = c("births", "deaths"), format = "csv") # One HMD-format .txt file per province download_ine_data("ine_data", stages = "life_tables", format = "txt") \end{lstlisting} \fn{"csv"} writes one combined file per data table (e.g.\ \fn{population.csv}, one row per province$\times$age$\times$year); \fn{"txt"} writes HMD-format files, one per province, matching the layout HMD-compatible tooling expects. \newpage \section{Part II --- Demographic Analysis and Visualization} \label{sec:part2} \subsection{Summary demographic indicators} The functions in this section all operate on System B's age/sex-disaggregated province data (or, for \fn{birth\_death\_ratio()}, on System A's totals --- noted explicitly below); \fn{get\_ine\_demog()}'s total-only indicators cannot be used for any indicator that requires an age breakdown. \subsubsection{Age dependency ratios: \fn{age\_dependency\_ratio()}} With $P_{\text{young}}$, $P_{\text{work}}$, $P_{\text{old}}$ the population aged $\leq 14$, $15$--$64$, and $\geq 65$ respectively (cut-offs configurable via \fn{young\_max}/\fn{old\_min}): \begin{align} \text{youth dependency ratio} &= \frac{P_{\text{young}}}{P_{\text{work}}}\times 100 \\ \text{old-age dependency ratio} &= \frac{P_{\text{old}}}{P_{\text{work}}}\times 100 \\ \text{total dependency ratio} &= \frac{P_{\text{young}}+P_{\text{old}}}{P_{\text{work}}}\times 100 \end{align} \begin{lstlisting} pop <- get_ine_population() adr <- age_dependency_ratio(pop$data) adr[adr$nuts3_code == "ES111", ] \end{lstlisting} \subsubsection{Aging index: \fn{aging\_index()}} \begin{equation} \text{aging index} = \frac{P_{\text{old}}}{P_{\text{young}}}\times 100 \end{equation} \begin{lstlisting} ai <- aging_index(pop$data) \end{lstlisting} \subsubsection{Sex ratio: \fn{sex\_ratio()}} \begin{equation} \text{sex ratio} = \frac{P_{\text{male}}}{P_{\text{female}}}\times 100 \end{equation} computed over all ages by default, or per single year of age with \fn{by\_age = TRUE}. \begin{lstlisting} sex_ratio(pop$data) sex_ratio(pop$data, by_age = TRUE) \end{lstlisting} \subsubsection{Crude birth and death rates} \begin{equation} \mathrm{CBR} = \frac{B}{P}\times 1000, \qquad \mathrm{CDR} = \frac{D}{P}\times 1000 \end{equation} where $B$, $D$, $P$ are total annual births, deaths, and population. \begin{lstlisting} births <- get_ine_births() deaths <- get_ine_deaths() cbr <- crude_birth_rate(births$data, pop$data) cdr <- crude_death_rate(deaths$data_provinces, pop$data) \end{lstlisting} \textbf{These are crude rates, not a total fertility rate (TFR).} TFR requires age-of-mother-specific birth counts; INE's underlying MNPN table as retrieved by \fn{get\_ine\_births()} carries no age-of-mother breakdown, and this package does not currently compute TFR. Do not substitute \fn{crude\_birth\_rate()} for TFR in fertility analysis. \subsubsection{General fertility rate: \fn{general\_fertility\_rate()}} GFR improves on CBR by dividing births by the female population actually at risk of childbearing (ages 15--49 by default), removing the distortion CBR suffers when two provinces have different age/sex structures despite similar underlying fertility behavior: \begin{equation} \mathrm{GFR} = \frac{B}{W_{15-49}} \times 1000 \end{equation} \begin{lstlisting} gfr <- general_fertility_rate(births$data, pop$data, age_min = 15, age_max = 49) \end{lstlisting} GFR is still not a TFR: it is a single aggregate rate over all reproductive-age women, not an age-specific fertility schedule. \subsubsection{Fertility schedule: ASFR, TFR, MAC, GRR, NRR} \label{sec:fertility-schedule} Where CBR and GFR (above) are single aggregate rates, the functions below use \fn{get\_ine\_births\_by\_age()}'s (Section~\ref{sec:births-by-age}) age-of-mother breakdown to build the full age-specific fertility schedule and the standard summary indicators derived from it. \paragraph{Age-specific fertility rate: \fn{age\_specific\_fertility\_rate()}.} \begin{equation} \mathrm{ASFR}(x) = \frac{B(x)}{W(x)} \times 1000, \qquad \mathrm{ASFR}_{\text{female}}(x) = \frac{B_{\text{female}}(x)}{W(x)} \times 1000 \end{equation} for single years of age $x = 15, \dots, 49$ by default, where $B(x)$ is births to mothers of age $x$ and $W(x)$ is the female population of age $x$. \fn{age\_specific\_fertility\_rate()} returns both the all-newborn rate $\mathrm{ASFR}(x)$ and the female-newborn-specific rate $\mathrm{ASFR}_{\text{female}}(x)$ (using INE's actual sex-of-newborn breakdown, not an assumed constant sex ratio at birth) --- the schedule the remaining functions below all consume. \begin{lstlisting} births_age <- get_ine_births_by_age() pop <- get_ine_population() asfr <- age_specific_fertility_rate(births_age$data, pop$data, age_min = 15, age_max = 49) \end{lstlisting} \subsubsection{Total fertility rate: \fn{total\_fertility\_rate()}} \begin{equation} \mathrm{TFR} = \sum_{x=15}^{49} \frac{\mathrm{ASFR}(x)}{1000} \end{equation} the true TFR that \fn{crude\_birth\_rate()} and \fn{general\_fertility\_rate()} explicitly note this package could not previously compute. \begin{lstlisting} tfr <- total_fertility_rate(asfr) \end{lstlisting} \subsubsection{Mean age at childbearing: \fn{mean\_age\_at\_childbearing()}} \begin{equation} \mathrm{MAC} = \frac{\sum_x x \cdot \mathrm{ASFR}(x)}{\sum_x \mathrm{ASFR}(x)} \end{equation} \begin{lstlisting} mac <- mean_age_at_childbearing(asfr) \end{lstlisting} \subsubsection{Gross and net reproduction rates} The gross reproduction rate (expected daughters per woman over her reproductive lifetime, ignoring mortality) sums the female-newborn-specific schedule directly, rather than approximating it by scaling TFR with an assumed constant sex ratio at birth: \begin{equation} \mathrm{GRR} = \sum_{x=15}^{49} \frac{\mathrm{ASFR}_{\text{female}}(x)}{1000} \end{equation} The net reproduction rate additionally weights that schedule by the mother's own survivorship to each reproductive age, using the female life table's $L_x$ (Section~\ref{sec:lifetables}) and radix $l_0$: \begin{equation} \mathrm{NRR} = \sum_{x=15}^{49} \frac{\mathrm{ASFR}_{\text{female}}(x)}{1000} \cdot \frac{L_x}{l_0} \end{equation} so $\mathrm{NRR} \leq \mathrm{GRR}$ always, with equality only in the (unrealistic) absence of any mortality before the end of the reproductive span. \begin{lstlisting} grr <- gross_reproduction_rate(asfr) nrr <- net_reproduction_rate(asfr, lt$fltper) # lt$fltper: female life table, build_life_tables() \end{lstlisting} As a live-verified example, Madrid 2023 gives $\mathrm{GRR} = 0.538$ and $\mathrm{NRR} = 0.535$ --- both well below the replacement threshold of 1, consistent with Spain's well-documented below-replacement fertility, with $\mathrm{NRR}$ correctly just under $\mathrm{GRR}$ given Spain's low mortality at reproductive ages. \subsubsection{Infant mortality rate: \fn{infant\_mortality\_rate()}} \begin{equation} \mathrm{IMR} = \frac{D_{0}}{B}\times 1000 \end{equation} where $D_0$ is deaths under age 1 in year $t$ and $B$ is live births in that \emph{same} year $t$ (the birth cohort those deaths are drawn from) --- a different denominator convention from CDR. IMR is closely related to, but not identical to, $q_0$ in the period life table: $q_0$ is derived from the exposure-based $m_0$ via the Andreev-Kingkade $a_0$, while IMR here uses the simpler, standard deaths/births ratio directly; the two will usually be close but need not match exactly. \begin{lstlisting} imr <- infant_mortality_rate(deaths$data_provinces, births$data) \end{lstlisting} \subsubsection{Rate of natural increase: \fn{rate\_of\_natural\_increase()}} \begin{equation} \mathrm{RNI} = \mathrm{CBR} - \mathrm{CDR} \end{equation} combining already-computed \fn{crude\_birth\_rate()} and \fn{crude\_death\_rate()} output. A positive value means the population is growing from natural change alone (ignoring migration); negative means the reverse --- the point in time and space where $\mathrm{CDR}$ overtakes $\mathrm{CBR}$ is often called the ``demographic crossover.'' \begin{lstlisting} rni <- rate_of_natural_increase(cbr, cdr) \end{lstlisting} \subsubsection{Birth-to-death ratio: \fn{birth\_death\_ratio()}} Unlike the indicators above, this needs no age breakdown, so it operates directly on \fn{get\_ine\_demog()}'s totals (System A): \begin{equation} \text{birth-death ratio} = \frac{B}{D} \end{equation} \begin{lstlisting} vitals <- get_ine_demog(indicator = c("births_total", "deaths_total"), year = 2023) ratio <- birth_death_ratio(vitals) \end{lstlisting} A value of 2 means 2 live births per death; 0.5 means 1 birth per 2 deaths. \subsubsection{Life expectancy: \fn{life\_expectancy\_summary()} and \fn{life\_expectancy()}} \fn{life\_expectancy\_summary()} extracts $e_0$ (life expectancy at birth) and $e_{65}$ (remaining life expectancy at 65) for every province and year from a \fn{build\_life\_tables()} output table: \begin{lstlisting} le <- life_expectancy_summary(lt$fltper, sex = "female") \end{lstlisting} \fn{life\_expectancy()} is a one-province convenience wrapper: it builds the life table for a single matched province directly from \fn{compute\_death\_rates()}'s \fn{mx\_1x1} output, without requiring \fn{build\_life\_tables()} to be run for every province first. \begin{lstlisting} life_expectancy(rates$mx_1x1, province = "A Coruna", sex = "female") \end{lstlisting} \subsection{Mortality decomposition: \fn{decompose\_life\_expectancy()}} \label{sec:decomposition} Life expectancy summaries answer \emph{what} $e_0$ is; mortality decomposition answers \emph{why} two life tables' $e_0$ values differ --- attributing the gap between two provinces (same year) or one province (two years) to age-specific contributions. Two independently-derived methods are implemented so results can be cross-checked against each other (Ponnapalli, 2005, found the two are not sensitive to which one is used): \begin{lstlisting} decompose_life_expectancy(lt1, lt2, method = c("arriaga", "pollard")) \end{lstlisting} \begin{description}[leftmargin=2.2cm,style=nextline] \item[lt1, lt2] Single life tables (one province-year-sex slice of \fn{build\_life\_tables()}'s \fn{fltper}/\fn{mltper}/\fn{bltper}, or \fn{build\_life\_table()}'s direct output), over the \emph{same} set of ages. \end{description} \paragraph{Arriaga's method (1984, the default).} An exact discrete decomposition. For a non-terminal age interval $[x, x+n)$, splitting into a direct effect (extra/fewer years lived \emph{within} the interval) and a combined indirect/interaction effect (extra/fewer survivors carrying forward into \emph{later} ages): \begin{align} \mathrm{DE}(x) &= \frac{l_x^{(1)}}{l_0^{(1)}} \left(\frac{{}_nL_x^{(2)}}{l_x^{(2)}} - \frac{{}_nL_x^{(1)}}{l_x^{(1)}}\right) \\ \mathrm{IE}(x) &= \frac{T_{x+n}^{(2)}}{l_0^{(1)}} \left(\frac{l_x^{(1)}}{l_x^{(2)}} - \frac{l_{x+n}^{(1)}}{l_{x+n}^{(2)}}\right) \\ \Delta(x) &= \mathrm{DE}(x) + \mathrm{IE}(x) \end{align} with superscripts $(1)$/$(2)$ denoting the reference/comparison life table. For the terminal open interval $\omega$, only a direct effect applies (since ${}_{\infty}L_\omega = T_\omega$): \begin{equation} \Delta(\omega) = \frac{l_\omega^{(1)}}{l_0^{(1)}} \left(\frac{T_\omega^{(2)}}{l_\omega^{(2)}} - \frac{T_\omega^{(1)}}{l_\omega^{(1)}}\right) \end{equation} and $\sum_x \Delta(x) = e_0^{(2)} - e_0^{(1)}$ exactly. \paragraph{Pollard's method (1988).} An exact \emph{continuous}-time decomposition; applied to a discrete single-year life table, its symmetric-average form (averaging the two individually-exact continuous identities Pollard derives, so the choice of reference table doesn't bias the result) is: \begin{equation} \Delta(x) = \tfrac{1}{2}\bigl(m_x^{(1)} - m_x^{(2)}\bigr) \left(\frac{l_x^{(2)}}{l_0^{(2)}} e_x^{(1)} + \frac{l_x^{(1)}}{l_0^{(1)}} e_x^{(2)}\right) \end{equation} Unlike Arriaga's, this is exact only in the continuous limit --- on a single-year life table it is a close but not identical approximation of $e_0^{(2)} - e_0^{(1)}$, so treat Arriaga's (the default) as the primary result and Pollard's as a cross-check. \begin{lstlisting} lt <- build_life_tables(rates$mx_1x1) lt_a <- lt$fltper[lt$fltper$province_name == "A Coruna" & lt$fltper$year == 2023, ] lt_b <- lt$fltper[lt$fltper$province_name == "Madrid" & lt$fltper$year == 2023, ] decomp <- decompose_life_expectancy(lt_a, lt_b, method = "arriaga") sum(decomp$contribution) # recovers the e0 gap exactly \end{lstlisting} As a live-verified example, A Coru\~{n}a vs.\ Madrid, 2023: $e_0$ gap $= 1.04$ years, with age 0 and age 58 the largest single-age contributors (Arriaga and Pollard correlate at $0.999$ on this comparison). For a temporal comparison --- A Coru\~{n}a, 2016 vs.\ 2023 --- the gap is $0.69$ years. \subsection{Mapping functions} \subsubsection{\fn{plot\_ine\_map()}} Plots the full map of Spain at province or municipality level, with one or more regions highlighted --- primarily a diagnostic tool for confirming which regions a regex \fn{region}/\fn{highlight} pattern actually matches. \begin{lstlisting} plot_ine_map(geo_level = "province", highlight = "^Cordoba$") \end{lstlisting} \subsubsection{\fn{map\_indicator()}} General-purpose choropleth: joins any tidy indicator tibble keyed by \fn{GEOID} onto province or municipality geometry, as either a binned discrete-legend map (default, matching the classic ``N regions above/below a threshold'' choropleth style) or a continuous gradient. \begin{lstlisting} map_indicator(df, value_col, geo_level = c("province", "municipality"), region = NULL, binned = TRUE, breaks = NULL, palette = if (binned) "PiYG" else "D", moveCAN = TRUE, can_gap_km = 60, legend_title = NULL, title = NULL) \end{lstlisting} \begin{lstlisting} pop_data <- get_ine_demog(indicator = "population_total", year = 2023) map_indicator(pop_data, "population_total", geo_level = "province") \end{lstlisting} \subsubsection{\fn{map\_life\_expectancy()}} A thin wrapper around \fn{map\_indicator()} that bridges System B's \fn{life\_expectancy\_summary()} output with System A's \fn{get\_ine\_geo()} geometry via the internal province lookup table. \begin{lstlisting} map_life_expectancy(le_df, year, sex = "both", moveCAN = TRUE, can_gap_km = 60, title = NULL) \end{lstlisting} \begin{lstlisting} map_life_expectancy(le, year = max(le$year), sex = "female") \end{lstlisting} \subsection{Population pyramids and trend charts} \subsubsection{\fn{plot\_population\_pyramid()}} Draws a mirrored population pyramid (female right, male left as negative counts) for one province/region and year. \begin{lstlisting} plot_population_pyramid(pop_df, year, region = NULL, title = NULL) \end{lstlisting} \begin{lstlisting} plot_population_pyramid(pop$data, year = max(pop$data$year), region = "A Coruna") \end{lstlisting} Note \fn{geom\_col(width = 1)} is used internally rather than \pkg{ggplot2}'s default \fn{width = 0.9}: at the default width, the systematic 10\% gap between each of the 101 single-year-age bars renders as visible thin white stripes across the pyramid. \subsubsection{\fn{plot\_demog\_trend()}} A generic time-series chart for any indicator tibble keyed by \fn{nuts3\_code}/\fn{province\_name}/\fn{year} --- works unmodified on the output of any indicator function above. \begin{lstlisting} plot_demog_trend(df, value_col, region = NULL, title = NULL, ylab = NULL) \end{lstlisting} \begin{lstlisting} plot_demog_trend(ai, "aging_index", region = "A Coruna") \end{lstlisting} \subsection{Lexis diagrams: \fn{plot\_lexis\_diagram()}} A Lexis diagram plots age (vertical axis) against calendar year (horizontal axis), shaded by mortality rate, with birth-cohort diagonals overlaid. Every point along one diagonal line represents the same birth cohort, aging one year for every calendar year that passes: \begin{equation} \text{age} = \text{year} - \text{cohort} \end{equation} \begin{lstlisting} plot_lexis_diagram(mx_df, province, sex = c("total", "female", "male"), log_scale = TRUE, rate_per = 1000, cohort_lines = TRUE, cohort_step = 10, title = NULL) \end{lstlisting} \begin{lstlisting} rates <- compute_death_rates(deaths$data_provinces, exposure$data) plot_lexis_diagram(rates$mx_1x1, province = "A Coruna") plot_lexis_diagram(rates$mx_1x1, province = "Madrid", cohort_lines = FALSE) \end{lstlisting} Mortality is color-scaled on a log10 axis by default (\fn{log\_scale = TRUE}), the standard convention for Lexis surfaces since $m_x$ spans several orders of magnitude across ages. Grey cells mean no rate could be computed for that age$\times$year cell --- almost always because the source population data has no single-year age breakdown for the oldest ages in the earliest years of the series (so exposure, the rate's denominator, is missing there), not a defect in the calculation; this typically appears as a rectangular block in the upper-left (oldest ages, earliest years) and clears once the source data's age detail becomes complete. When present, grey is given its own ``Missing data'' legend key via a small dummy point layer, since a continuous color scale cannot show an \fn{na.value} swatch on its own. \subsection{Worked examples} The nine scripts summarized below live in \fn{examples/} and are fully runnable end to end (each retrieves live data with the package's own functions before plotting). They double as integration tests for the indicator and charting functions documented above. \subsubsection{Birth-to-death ratio map (\fn{birth\_death\_ratio\_map.R})} Replicates a reference choropleth style (binned, pink-to-green diverging scale) at the finest granularity Spain's public birth/death data actually supports: province (municipality-level birth/death tables do not exist in INE's MNPN/MNPD operations). \begin{lstlisting} library(inedemogR) vitals <- get_ine_demog(indicator = c("births_total", "deaths_total"), year = 2023) ratio <- birth_death_ratio(vitals) map_indicator( ratio, value_col = "birth_death_ratio", geo_level = "province", breaks = c(0.25, 0.5, 0.75, 1, 1.5, 2, 3), palette = "PiYG", legend_title = "Births per death", title = "Birth-to-death ratio by province, Spain (2023)" ) \end{lstlisting} \subsubsection{Andalusian population pyramids (\fn{andalusia\_pyramids.R})} Builds an eight-panel \fn{facet\_wrap()} grid of population pyramids for Andalusia's eight provinces (all sharing \fn{nuts2\_code == "ES61"}), directly in \pkg{ggplot2} syntax rather than calling \fn{plot\_population\_pyramid()} eight times. \begin{lstlisting} library(inedemogR); library(dplyr); library(tidyr); library(ggplot2) pop <- get_ine_population() latest_year <- max(pop$data$year) andalusia <- pop$data |> filter(nuts2_code == "ES61", year == latest_year) pyramid_data <- andalusia |> select(province_name, age, female, male) |> pivot_longer(c(female, male), names_to = "sex", values_to = "count") |> mutate(count = if_else(sex == "male", -count, count)) ggplot(pyramid_data, aes(x = age, y = count, fill = sex)) + geom_col(width = 1) + coord_flip() + facet_wrap(~province_name, nrow = 2, ncol = 4, scales = "free_x") + scale_y_continuous(labels = abs) + scale_fill_manual(values = c(female = "#D4667A", male = "#4477AA")) + labs(x = "Age", y = "Population", fill = "Sex", title = paste("Population pyramids, Andalusian provinces,", latest_year)) + theme_minimal() \end{lstlisting} \subsubsection{Madrid/Barcelona Lexis diagrams (\fn{madrid\_barcelona\_lexis.R})} \fn{plot\_lexis\_diagram()} handles one province per call; a two-panel comparison is built directly with the pipeline functions plus \fn{facet\_wrap()}. \begin{lstlisting} library(inedemogR); library(dplyr); library(ggplot2) pop <- get_ine_population() deaths <- get_ine_deaths() exposure <- compute_exposure(pop$data, deaths$data_provinces) rates <- compute_death_rates(deaths$data_provinces, exposure$data) mx <- rates$mx_1x1 |> filter(province_name %in% c("Madrid", "Barcelona")) |> mutate(mx_per_1000 = mx_total * 1000) ggplot(mx, aes(x = year, y = age, fill = mx_per_1000)) + geom_tile() + scale_fill_viridis_c( name = "Mortality rate\n(per 1,000, log scale)", trans = "log10", na.value = "grey50", labels = scales::label_number(accuracy = 0.1) ) + facet_wrap(~province_name, nrow = 1, ncol = 2) + labs(x = "Year", y = "Age", title = "Lexis diagrams: Madrid and Barcelona") + theme_minimal() \end{lstlisting} \subsubsection{Life expectancy comparison (\fn{life\_expectancy\_comparison.R})} Combines a multi-province time series with a national choropleth for one year, both from the same life-table pipeline run once for every province. \begin{lstlisting} library(inedemogR); library(dplyr); library(ggplot2) pop <- get_ine_population() deaths <- get_ine_deaths() exposure <- compute_exposure(pop$data, deaths$data_provinces) rates <- compute_death_rates(deaths$data_provinces, exposure$data) lt <- build_life_tables(rates$mx_1x1) le <- life_expectancy_summary(lt$fltper, sex = "female") provinces_to_compare <- c("Madrid", "Barcelona", "Sevilla", "A Coruna") ggplot(le |> filter(province_name %in% provinces_to_compare), aes(x = year, y = e0, color = province_name)) + geom_line(linewidth = 0.8) + geom_point(size = 1.5) + labs(x = "Year", y = "Life expectancy at birth (e0)", color = "Province") + theme_minimal() map_life_expectancy(le, year = max(le$year), sex = "female", title = "Life expectancy at birth by province, Spain") \end{lstlisting} \subsubsection{Natural increase analysis (\fn{natural\_increase\_analysis.R})} Visualizes the ``demographic crossover'' (CDR overtaking CBR) and the resulting rate of natural increase. \begin{lstlisting} library(inedemogR); library(dplyr); library(ggplot2) births <- get_ine_births(); deaths <- get_ine_deaths(); pop <- get_ine_population() cbr <- crude_birth_rate(births$data, pop$data) cdr <- crude_death_rate(deaths$data_provinces, pop$data) rni <- rate_of_natural_increase(cbr, cdr) provinces_to_compare <- c("Madrid", "Barcelona", "Sevilla", "A Coruna") ggplot(rni |> filter(province_name %in% provinces_to_compare), aes(x = year, y = rni, color = province_name)) + geom_hline(yintercept = 0, linetype = "dashed", color = "grey50") + geom_line(linewidth = 0.8) + geom_point(size = 1.5) + labs(x = "Year", y = "Rate of natural increase (per 1,000)", color = "Province") + theme_minimal() \end{lstlisting} \subsubsection{Fertility and infant mortality analysis (\fn{fertility\_infant\_mortality\_analysis.R})} \begin{lstlisting} library(inedemogR); library(dplyr); library(ggplot2) births <- get_ine_births(); deaths <- get_ine_deaths(); pop <- get_ine_population() gfr <- general_fertility_rate(births$data, pop$data) imr <- infant_mortality_rate(deaths$data_provinces, births$data) provinces_to_compare <- c("Madrid", "Barcelona", "Sevilla", "A Coruna") ggplot(gfr |> filter(province_name %in% provinces_to_compare), aes(x = year, y = gfr, color = province_name)) + geom_line(linewidth = 0.8) + geom_point(size = 1.5) + labs(x = "Year", y = "General fertility rate (per 1,000 women 15-49)", color = "Province") + theme_minimal() \end{lstlisting} \subsubsection{Fertility indicators analysis (\fn{fertility\_indicators\_analysis.R})} Exercises the full ASFR-based fertility schedule (Section~\ref{sec:fertility-schedule}): the age schedule itself, TFR/MAC trends, and the GRR-vs-NRR mortality-adjustment gap, for a comparison set of provinces. \begin{lstlisting} library(inedemogR); library(dplyr); library(ggplot2) births_age <- get_ine_births_by_age() pop <- get_ine_population() deaths <- get_ine_deaths() asfr <- age_specific_fertility_rate(births_age$data, pop$data) tfr <- total_fertility_rate(asfr) mac <- mean_age_at_childbearing(asfr) grr <- gross_reproduction_rate(asfr) exposure <- compute_exposure(pop$data, deaths$data_provinces) rates <- compute_death_rates(deaths$data_provinces, exposure$data) lt <- build_life_tables(rates$mx_1x1) nrr <- net_reproduction_rate(asfr, lt$fltper) latest_year <- max(asfr$year) provinces_to_compare <- c("Madrid", "Barcelona", "Sevilla", "A Coruna") ggplot(asfr |> filter(province_name %in% provinces_to_compare, year == latest_year), aes(x = age, y = asfr, color = province_name)) + geom_line(linewidth = 0.8) + geom_point(size = 1.2) + labs(x = "Age of mother", y = "ASFR (births per 1,000 women of that age)", color = "Province", title = paste("ASFR schedule,", latest_year)) + theme_minimal() ggplot(tfr |> filter(province_name %in% provinces_to_compare), aes(x = year, y = tfr, color = province_name)) + geom_line(linewidth = 0.8) + geom_point(size = 1.5) + geom_hline(yintercept = 2.1, linetype = "dashed", color = "grey40") + labs(x = "Year", y = "Total fertility rate", color = "Province", subtitle = "Dashed line: replacement-level TFR (2.1)") + theme_minimal() \end{lstlisting} Live-verified results: TFR ranges roughly $1.0$--$1.3$ across the four provinces (Madrid highest, A Coru\~{n}a lowest), well below replacement, and Madrid's mean age at childbearing comes to $32.8$ years --- consistent with Spain's known late-childbearing pattern. \subsubsection{Mortality decomposition analysis (\fn{mortality\_decomposition\_analysis.R})} Exercises \fn{decompose\_life\_expectancy()} (Section~\ref{sec:decomposition}) on both a cross-sectional gap (two provinces, one year) and a temporal gap (one province, two years), comparing Arriaga's and Pollard's methods. \begin{lstlisting} library(inedemogR); library(dplyr); library(ggplot2) pop <- get_ine_population(); deaths <- get_ine_deaths() exposure <- compute_exposure(pop$data, deaths$data_provinces) rates <- compute_death_rates(deaths$data_provinces, exposure$data) lt <- build_life_tables(rates$mx_1x1) latest_year <- max(lt$fltper$year); earliest_year <- min(lt$fltper$year) lt_a <- lt$fltper |> filter(province_name == "A Coruna", year == latest_year) lt_b <- lt$fltper |> filter(province_name == "Madrid", year == latest_year) decomp <- decompose_life_expectancy(lt_a, lt_b, method = "arriaga") ggplot(decomp, aes(x = age, y = contribution, fill = contribution > 0)) + geom_col(width = 1) + scale_fill_manual(values = c(`TRUE` = "steelblue", `FALSE` = "firebrick"), guide = "none") + labs(x = "Age", y = "Contribution to e0 gap (years)", title = "Madrid vs. A Coruna life expectancy gap by age") + theme_minimal() \end{lstlisting} Live-verified: A Coru\~{n}a 2024 $e_0 = 86.97$ (matching the value confirmed in Section~\ref{sec:lifetables} exactly), Madrid 2024 $e_0 = 87.83$ --- gap $0.86$ years; the temporal comparison (A Coru\~{n}a, 1996--2024) gives a gap of $5.23$ years, consistent with the secular mortality improvement documented throughout this tutorial. \subsubsection{Abridged life tables analysis (\fn{abridged\_life\_tables\_analysis.R})} Cross-checks \fn{build\_abridged\_life\_tables()} (Section~\ref{sec:abridged}) against the exact single-year \fn{build\_life\_tables()} on the same underlying rates, then charts the abridged survivorship curve and mortality schedule. \begin{lstlisting} library(inedemogR); library(dplyr); library(ggplot2) pop <- get_ine_population(); deaths <- get_ine_deaths() exposure <- compute_exposure(pop$data, deaths$data_provinces) rates <- compute_death_rates(deaths$data_provinces, exposure$data) lt_1x1 <- build_life_tables(rates$mx_1x1) lt_abridged <- build_abridged_life_tables(rates$mx_1x1, rates$mx_5x1) latest_year <- max(lt_abridged$fltper$year) lt_1x1_one <- lt_1x1$fltper |> filter(province_name == "A Coruna", year == latest_year) lt_abridged_one <- lt_abridged$fltper |> filter(province_name == "A Coruna", year == latest_year) survivorship <- bind_rows( lt_1x1_one |> transmute(age = age, lx = lx, table = "Single-year (1x1)"), lt_abridged_one |> transmute(age = age_start, lx = lx, table = "Abridged (5x1)") ) ggplot(survivorship, aes(x = age, y = lx, color = table)) + geom_line(linewidth = 0.8) + labs(x = "Age", y = "Survivors (out of 100,000 births)", color = NULL) + theme_minimal() \end{lstlisting} Live-verified: A Coru\~{n}a 2024 abridged $e_0 = 87.03$ against the exact single-year $86.97$ --- a $0.06$-year difference, confirming the abridged construction recovers the single-year result closely. \newpage \appendix \section*{Appendix} \addcontentsline{toc}{section}{Appendix} \subsection*{A.1 Function reference} \begin{longtable}{@{}p{5.3cm}p{9.5cm}@{}} \toprule Function & Purpose \\ \midrule \endhead \fn{get\_ine\_demog()} & System A: totals for one or more indicators, any geography \\ \fn{get\_ine\_geo()} & Province/municipality boundary geometries \\ \fn{list\_ine\_indicators()} & The indicator registry \\ \fn{update\_ine\_data()} & Refresh cached System A indicators if INE published newer data \\ \fn{plot\_ine\_map()} & Highlight-region diagnostic map \\ \fn{map\_indicator()} & General choropleth for any \fn{GEOID}-keyed tibble \\ \fn{get\_ine\_population()} & Province population by single-year age/sex \\ \fn{get\_ine\_births()} & Province births by sex \\ \fn{get\_ine\_births\_by\_age()} & Province births by single year of age of mother and sex of newborn \\ \fn{get\_ine\_deaths()} & Province deaths by sex, age-specific and age-less \\ \fn{compute\_exposure()} & Exposure-to-risk $E(x,t)$ \\ \fn{compute\_death\_rates()} & Central death rates: 1x1, 5x1, ASDR \\ \fn{andreev\_kingkade\_a0()} & Age-0 life-table parameter $a_0$ \\ \fn{build\_life\_table()} / \fn{build\_life\_tables()} & Period life table, one sex/one province, or all \\ \fn{build\_abridged\_life\_table()} / \fn{build\_abridged\_life\_tables()} & Abridged (5-year age group) period life table, one sex/one province, or all \\ \fn{download\_ine\_data()} & Export the whole pipeline to CSV/HMD-txt files \\ \fn{validate\_population()}, \fn{validate\_births()}, \fn{validate\_births\_by\_age()}, \fn{validate\_deaths()}, \fn{validate\_deaths\_age()}, \fn{validate\_exposure()}, \fn{validate\_death\_rates()}, \fn{validate\_life\_table()}, \fn{validate\_abridged\_life\_table()} & Quality-control checks for each pipeline stage \\ \fn{age\_dependency\_ratio()} & Youth/old-age/total dependency ratios \\ \fn{aging\_index()} & Elderly per 100 children \\ \fn{sex\_ratio()} & Males per 100 females \\ \fn{crude\_birth\_rate()} / \fn{crude\_death\_rate()} & CBR / CDR \\ \fn{general\_fertility\_rate()} & GFR (births per 1,000 women 15-49) \\ \fn{age\_specific\_fertility\_rate()} & ASFR/ASFR-female schedule by single year of age \\ \fn{total\_fertility\_rate()} & TFR (children per woman) \\ \fn{mean\_age\_at\_childbearing()} & MAC \\ \fn{gross\_reproduction\_rate()} & GRR (daughters per woman, ignoring mortality) \\ \fn{net\_reproduction\_rate()} & NRR (daughters per woman, mortality-adjusted) \\ \fn{infant\_mortality\_rate()} & IMR (infant deaths per 1,000 live births) \\ \fn{rate\_of\_natural\_increase()} & RNI = CBR $-$ CDR \\ \fn{birth\_death\_ratio()} & Births per death, System A totals \\ \fn{life\_expectancy\_summary()} & $e_0$/$e_{65}$ from a full life-table output \\ \fn{life\_expectancy()} & $e_0$/$e_{65}$, one-province convenience wrapper \\ \fn{decompose\_life\_expectancy()} & Age-specific decomposition of an $e_0$ gap (Arriaga/Pollard) \\ \fn{plot\_population\_pyramid()} & Mirrored age/sex pyramid \\ \fn{plot\_demog\_trend()} & Generic indicator time series \\ \fn{plot\_lexis\_diagram()} & Age $\times$ year mortality surface with cohort diagonals \\ \fn{map\_life\_expectancy()} & Life-expectancy choropleth (bridges System B to System A geometry) \\ \bottomrule \end{longtable} \subsection*{A.2 INE data sources used} \begin{longtable}{@{}p{3.3cm}p{2cm}p{9.5cm}@{}} \toprule Series & Table ID & Operation \\ \midrule \endhead Population (municipality) & 29005 & Cifras oficiales del padr\'{o}n (DPOP) \\ Population (province, by age/sex) & 56945 & Padr\'{o}n Municipal Continuo (ECP) \\ Births (province) & 6506 & Movimiento Natural de la Poblaci\'{o}n --- Nacimientos (MNPN) \\ Births (province, age of mother) & 6508 & Movimiento Natural de la Poblaci\'{o}n --- Nacimientos (MNPN) \\ Deaths (province, age-less) & 6545 & Movimiento Natural de la Poblaci\'{o}n --- Defunciones (MNPD) \\ Deaths (province, age-specific) & 6547 & Movimiento Natural de la Poblaci\'{o}n --- Defunciones (MNPD) \\ \bottomrule \end{longtable} \subsection*{A.3 Methodology references} \begin{itemize} \item Human Mortality Database. \emph{Methods Protocol for the Human Mortality Database}, Version 6. University of California, Berkeley, and Max Planck Institute for Demographic Research. \item Andreev, E.\ M., \& Kingkade, W.\ W.\ (2015). Average age at death in infancy and infant mortality level: Reconsidering the Coale-Demeny formulas at current levels of low mortality. \emph{Demographic Research}, 33, 363--390. \item Eurostat. \emph{Revision of the European Standard Population} (2013 edition). \item Preston, S.\ H., Heuveline, P., \& Guillot, M.\ (2001). \emph{Demography: Measuring and Modeling Population Processes}. Blackwell Publishing. (Abridged life-table construction, Ch.\ 3.) \item Arriaga, E.\ E.\ (1984). Measuring and explaining the change in life expectancies. \emph{Demography}, 21(1), 83--96. \item Pollard, J.\ H.\ (1988). On the decomposition of changes in expectation of life and differentials in life expectancy. \emph{Demography}, 25(2), 265--276. \item Instituto Nacional de Estad\'{i}stica (INE). \emph{Fen\'{o}menos demogr\'{a}ficos}. \url{https://www.ine.es}. \end{itemize} \end{document}