%\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$. HMD Methods Protocol V6 (Eq.~57, uniform
distribution of births within cohorts) gives

\begin{equation}
E(x,t) \;=\; \tfrac{1}{2}\bigl[P(x,t) + P(x,t+1)\bigr]
\;+\; \tfrac{1}{6}\bigl[D^{L}(x,t) - D^{U}(x,t)\bigr]
\end{equation}

where $P(x,t)$ is the January~1 population aged $x$ in year $t$;
$D^{L}(x,t)$ counts deaths in the lower Lexis triangle (the cohort born in
$t-x$, which reaches age $x$ during year $t$) and $D^{U}(x,t)$ deaths in the
upper triangle (the cohort born in $t-x-1$, already aged $x$ on 1~January).
If \fn{deaths} carries a \fn{cohort} column the triangles are taken from it.
INE's death tables are not classified by cohort, so by default deaths are
split evenly between the triangles, the correction vanishes, and exposure is
the mean of the two January-1 stocks. This is an \emph{adaptation} of the
protocol, which splits square-cell deaths by regression (Appendix~A) and
corrects for the monthly distribution of births (Appendix~E); the
sensitivity of $e_0$ and $e_{65}$ to the split is reported in the package's
validation study (\texttt{reproducibility/}).

\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.
If $P(x,t+1)$ is missing for a year that is kept, exposure falls back to
$P(x,t)$ and the row is flagged \fn{is\_boundary\_year}.

\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,
any remaining negative or missing exposure, and exposure exceeding twice the
source population.

\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}

\paragraph{Missing data.} Three situations are kept apart:
\begin{itemize}
  \item \emph{Structural zeros.} INE's age-specific deaths table omits cells
  with no deaths, so an age missing from a province-year that is otherwise
  present is a true zero.
  \item \emph{Unavailable or suppressed counts.} An explicit \fn{NA} death
  count (INE suppresses some small sex-specific cells; those recoverable as
  Total minus the other sex are recovered at retrieval) stays \fn{NA} and
  propagates to the rate, the grouped rate, and the ASDR, and makes
  \fn{validate\_death\_rates()} fail.
  \item \emph{Missing province-years.} A province-year with exposure but no
  death records at all is excluded with a warning; it is never treated as
  zero mortality.
\end{itemize}

A pooled 5-year-age-group ($n=5$) rate is also computed from the same counts:

\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-standardised death rate per 100{,}000 using the European Standard
Population 2013 (ESP2013):

\begin{equation}
\mathrm{ASDR}(t) \;=\; \sum_{x=0}^{100+} m(x,t)\, w(x), \qquad \sum_x w(x) = 100{,}000.
\end{equation}

ESP2013 gives weights to five-year groups up to an open group 95+ (weight
200). Each group's weight is spread evenly over its single ages; the 95+
weight is spread over the six cells $95,\dots,99,100+$. The ASDR is \fn{NA}
unless every age has a finite rate.

\begin{lstlisting}
compute_death_rates(deaths, exposure)
\end{lstlisting}

Returns \fn{list(mx\_1x1, mx\_5x1, asdr, qc)}; \fn{mx\_1x1} and
\fn{mx\_5x1} carry the deaths (\fn{d\_*}) and exposures (\fn{e\_*}) beside
the rates, as the old-age model below needs them.
\fn{validate\_death\_rates()} fails on negative, missing or infinite rates,
on rates $\geq 1$ below age 90, and on incomplete age coverage.

\begin{lstlisting}
rates <- compute_death_rates(deaths$data_provinces, exposure$data)
rates$mx_1x1   # age x year x province mortality schedule
rates$asdr     # age-standardised rate per 100,000, one row per province-year
\end{lstlisting}

\subsubsection{Period life tables: \fn{build\_life\_tables()}}
\label{sec:lifetables}

\fn{build\_life\_tables()} constructs single-year period life tables
(ages 0 to 110+) for every province and year in \fn{mx\_1x1}, following
HMD Methods Protocol V6, Section~7.1.

\paragraph{Kannisto smoothing of old-age mortality.} Observed rates are
replaced from the age $Y$ upward, where $Y$ is the lowest age from 80 with
at most 100 female or 100 male deaths, constrained to $80 \le Y \le 95$
(the same $Y$ for both sexes). The fitted hazard is the Kannisto logistic
model with asymptote one,

\begin{equation}
\mu_x(a,b) = \frac{a\,e^{b(x-80)}}{1 + a\,e^{b(x-80)}}, \qquad a, b \ge 0,
\end{equation}

estimated by Poisson maximum likelihood, $D_x \sim \mathrm{Poisson}(E_x\,
\mu_{x+1/2})$, maximising
$\sum_x \bigl[D_x \log \mu_{x+1/2} - E_x\, \mu_{x+1/2}\bigr]$ over ages
80--99 (L-BFGS-B, analytic gradient, grid-search start). HMD fits ages
80--110+; INE's single-year detail ends at an open 100+ group, so the
package fits 80--99. The constraint $b \ge 0$ means smoothed rates cannot
fall with age. The smoothed rate at age $x$ is $\hat\mu_{x+1/2}$.
For the combined-sex table, rates from $Y$ are the sex-specific smoothed
rates weighted by a smoothed female share of exposure (weighted least
squares of $\mathrm{logit}\,\pi^F_x$ on a quadratic in age, HMD Eqs.~69--73).

\paragraph{Andreev--Kingkade $a_0$.} The mean time lived in the first year
by those who die before age 1 uses the period formulas of Andreev and
Kingkade (2015) as tabulated in HMD Methods Protocol V6, Table~3:

\begin{equation}
a_0 =
\begin{cases}
0.14929 - 1.99545\,m_0 & \text{male}, \ m_0 < 0.02300 \\
0.02832 + 3.26021\,m_0 & \text{male}, \ 0.02300 \le m_0 < 0.08307 \\
0.29915 & \text{male}, \ m_0 \ge 0.08307 \\
0.14903 - 2.05527\,m_0 & \text{female}, \ m_0 < 0.01724 \\
0.04667 + 3.88089\,m_0 & \text{female}, \ 0.01724 \le m_0 < 0.06891 \\
0.31411 & \text{female}, \ m_0 \ge 0.06891
\end{cases}
\end{equation}

For the combined-sex table, $a_0$ is the average of the female and male
values weighted by their age-0 deaths (HMD Eq.~77). Every other closed age
has $a_x = 0.5$.

\paragraph{Life-table recursion.} Given $m_x$ and $a_x$ for $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, \qquad d_x = l_x\, q_x, \qquad 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, \qquad e_x = T_x / l_x
\end{align}

so that $m_x = d_x / L_x$ holds at every age, including the open interval
($a_\omega = 1/m_\omega$).

\paragraph{Failure behaviour and coverage.} A province-year is tabulated only
if exposure is available at every single age 0--100 and rates are finite
below $Y$; tables that do not qualify are withheld and listed, with the
reason, in \fn{\$failed}. INE's provincial population by single year of age
is top-coded at 85+ before 2002, so the supported period is 2002 onwards.
Zero exposure at ages $\ge Y$ (for example, very old men in Ceuta and
Melilla) is allowed, because those rates come from the fitted model.
\fn{validate\_life\_table()} fails on missing or non-contiguous ages,
non-finite values, increasing $l_x$, non-positive $L_x$ or $e_x$,
$q_x \notin [0,1]$, a broken $m_x = d_x/L_x$ identity, or
$e_0 \notin (15, 100]$. The arguments \fn{kannisto\_age} (a fixed $Y$) and
\fn{fit\_min} (first fitting age) support sensitivity analysis.

\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,
failed)}: 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}; a \fn{qc}
tibble (one row per table, with the replacement age $Y$); and the withheld
tables.

\begin{lstlisting}
lt <- build_life_tables(rates$mx_1x1)
lt$fltper[lt$fltper$age == 0, ]   # e0 by province and year, female
\end{lstlisting}

With the frozen inputs of the package's \texttt{reproducibility/} folder,
this gives for A Coru\~{n}a, 2024, $e_0 = 86.95$ years (female) and $81.21$
years (male).

\subsubsection{Abridged life tables: \fn{build\_abridged\_life\_tables()}}
\label{sec:abridged}

\fn{build\_abridged\_life\_tables()} builds abridged period life tables
(\fn{"00-04"}, \fn{"05-09"}, \dots, \fn{"95-99"}, \fn{"100+"}) directly from
the observed grouped rates \fn{mx\_5x1}. It is an independent estimate, not
an abridgement of the single-year table: it uses no Kannisto smoothing and
closes at 100+ rather than 110+.

The youngest group is the one exception to using \fn{mx\_5x1} alone: a
single-year sub-table for ages 0--4 (Andreev--Kingkade $a_0$, $a_x = 0.5$
at ages 1--4) is built from \fn{mx\_1x1} and aggregated, which avoids an
external ${}_5a_0$ approximation for the sharp infant gradient.

For every other closed interval $[x, x+n)$ ($n=5$) the package uses
Greville's ${}_na_x$ with the matching Chiang conversion (Preston, Heuveline
\& Guillot, 2001, Box~3.1):

\begin{align}
{}_na_x &= \frac{n}{2} - \frac{n^2}{12}\Bigl({}_nm_x - k_x\Bigr), \qquad
k_x = \frac{1}{2n}\log\frac{{}_nm_{x+n}}{{}_nm_{x-n}} \\
{}_nq_x &= \frac{n\,{}_nm_x}{1 + (n - {}_na_x)\,{}_nm_x}, \qquad
{}_nL_x = n\, l_{x+n} + {}_na_x\, {}_nd_x
\end{align}

This pair satisfies ${}_nm_x = {}_nd_x / {}_nL_x$ exactly. Where $k_x$ is
undefined (a zero rate in a neighbouring group) or ${}_nq_x$ would exceed 1,
the group falls back to a constant hazard,
${}_nq_x = 1 - e^{-n\,{}_nm_x}$ and ${}_nL_x = {}_nd_x/{}_nm_x$, which also
preserves the identity. The open interval has $q = 1$ and $L = l/m$.

\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()} also
checks the identity ${}_nm_x = {}_nd_x/{}_nL_x$ in every closed group.

\begin{lstlisting}
alt <- build_abridged_life_tables(rates$mx_1x1, rates$mx_5x1)
alt$fltper[alt$fltper$age_group == "00-04", ]
\end{lstlisting}

For A Coru\~{n}a, 2024, the abridged female $e_0$ is $86.95$ years, within
$0.01$ years of the single-year table. Because the two tables share inputs
and part of their method, this agreement is a consistency check, not an
external validation; the comparison with INE's published tables is reported
in \texttt{reproducibility/}.

\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()}.}
The package follows INE's methodology for its Indicadores Demogr\'aficos
B\'asicos:

\begin{equation}
\mathrm{ASFR}(x) = \frac{B(x)}{\bar W(x)}, \qquad
\mathrm{ASFR}_{\text{female}}(x) = \frac{B_{\text{female}}(x)}{\bar W(x)}, \qquad
\bar W(x) = \tfrac{1}{2}\bigl[W(x,t) + W(x,t+1)\bigr]
\end{equation}

for single years of age $x = 15, \dots, 49$, where $B(x)$ is births in year
$t$ to mothers of completed age $x$ and $\bar W(x)$ is the mean annual
female population (the mean of the January-1 stocks of $t$ and $t+1$; if
the $t+1$ stock is not available the $t$ stock is used and the row is
flagged \fn{is\_boundary\_year}). Births to mothers younger than 15 are added
to age 15 and births to mothers aged 50 or over to age 49, so no births are
lost. Rates are births \emph{per woman}; the column \fn{asfr\_per\_1000}
gives the same schedule per 1{,}000 women for display.
$\mathrm{ASFR}_{\text{female}}$ uses INE's sex-of-newborn breakdown rather
than an assumed sex ratio at birth.

\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} \mathrm{ASFR}(x) \quad \text{(children per woman)}
\end{equation}

\begin{lstlisting}
tfr <- total_fertility_rate(asfr)
\end{lstlisting}

\subsubsection{Mean age at childbearing: \fn{mean\_age\_at\_childbearing()}}

Completed age $x$ covers the interval $[x, x+1)$, so its mid-point is used:

\begin{equation}
\mathrm{MAC} = \frac{\sum_x (x + 0.5)\, \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}

\begin{equation}
\mathrm{GRR} = \sum_{x=15}^{49} \mathrm{ASFR}_{\text{female}}(x), \qquad
\mathrm{NRR} = \sum_{x=15}^{49} \mathrm{ASFR}_{\text{female}}(x)\, \frac{L_x}{l_0}
\end{equation}

using the female period life table of the same province and year
(Section~\ref{sec:lifetables}), so $\mathrm{NRR} \leq \mathrm{GRR}$.

\begin{lstlisting}
grr <- gross_reproduction_rate(asfr)
nrr <- net_reproduction_rate(asfr, lt$fltper)   # lt$fltper: female life table, build_life_tables()
\end{lstlisting}

With the frozen inputs, Madrid 2023 gives $\mathrm{GRR} = 0.533$ and
$\mathrm{NRR} = 0.529$.

\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. The package uses the symmetric average of Pollard's two exact
forms,
$e_0^{(2)} - e_0^{(1)} = \int \bigl(\mu^{(1)}(x) - \mu^{(2)}(x)\bigr)\,
\tfrac12\bigl[\tfrac{l^{(2)}(x)}{l_0^{(2)}} e^{(1)}(x) + \tfrac{l^{(1)}(x)}{l_0^{(1)}} e^{(2)}(x)\bigr]\,dx$,
and integrates each closed single-year interval by the midpoint rule,

\begin{equation}
\Delta(x) = \tfrac{1}{2}\bigl(m_x^{(1)} - m_x^{(2)}\bigr)
  \left(\frac{L_x^{(2)}}{l_0^{(2)}}\,\bar e_x^{(1)} + \frac{L_x^{(1)}}{l_0^{(1)}}\,\bar e_x^{(2)}\right),
  \qquad \bar e_x = \tfrac12\,(e_x + e_{x+1}),
\end{equation}

and the open interval $\omega$ exactly under its constant hazard
($l(\omega+t) = l_\omega e^{-\mu t}$, $e(\omega+t) = 1/\mu$):

\begin{equation}
\Delta(\omega) = \bigl(m_\omega^{(1)} - m_\omega^{(2)}\bigr)\, e_\omega^{(1)} e_\omega^{(2)}
  \cdot \tfrac12\left(\frac{l_\omega^{(1)}}{l_0^{(1)}} + \frac{l_\omega^{(2)}}{l_0^{(2)}}\right).
\end{equation}

The sum differs from $e_0^{(2)} - e_0^{(1)}$ by a small discretisation
residual, returned with the result:

\begin{lstlisting}
lt <- build_life_tables(rates$mx_1x1)
female <- lt$fltper   # female tables
lt_a <- female[female$province_name == "A Coruna" & female$year == 2023, ]
lt_b <- female[female$province_name == "Madrid" & female$year == 2023, ]

decomp <- decompose_life_expectancy(lt_a, lt_b, method = "arriaga")
sum(decomp$contribution)          # recovers the e0 gap exactly
attr(decomp, "e0_difference")     # e0(Madrid) - e0(A Coruna)
attr(decompose_life_expectancy(lt_a, lt_b, "pollard"), "residual")
\end{lstlisting}

With the frozen inputs, the female $e_0$ gap between A Coru\~{n}a and
Madrid in 2023 is $1.04$ years, with ages 0, 58 and 90 the largest
single-age contributors; Pollard's contributions sum to the gap within
$0.0003$ years. For a temporal comparison, A Coru\~{n}a 2016 versus 2023, the
gain is $0.76$ 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 = 2024, 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 = 2024, 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 <- 2024
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 = 2024, 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 <- 2024
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_per_1000, 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 <- 2024; earliest_year <- 2006   # fixed years: results do not change with new INE releases
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}

With the frozen inputs: A Coru\~{n}a 2024 female $e_0 = 86.95$, Madrid 2024
$e_0 = 87.86$, a gap of $0.91$ years; the temporal comparison (A Coru\~{n}a,
2006--2024) gives a gain of $2.85$ years.

\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 <- 2024
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}

With the frozen inputs, A Coru\~{n}a 2024 abridged female $e_0 = 86.95$,
within $0.01$ years of the single-year table. The two tables share their
inputs, so this is a consistency check rather than external validation.

\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}