Package 'bage'

bage: Bayesian Estimation and Forecasting of Age-Specific Rates

Description

Modeling of rates, probabilities, and other values, typically disaggregated by age. Estimation is done using TMB, which makes it fast and scalable.

Example workflow

1. Specify model using mod_pois()

2. Fit model using fit()

3. Extract results using augment()

4. Check model using replicate_data()

Functions

Specify model

• mod_pois() Specify a Poisson model

• mod_binom() Specify a binomial model

• mod_norm() Specify a normal model

• set_prior() Specify prior for main effect or interaction

• priors Overview of priors for main effects or interactions

• set_disp() Specify prior for dispersion/variance

• set_covariates() Add covariates to model

• datamods Overview of data models (measurement error models)

• confidential Overview of confidentialization models

Fit model

• fit() Derive posterior distribution

Extract results

• augment() Original data, plus observation-level estimates

• components() Hyper-parameters

• dispersion() Dispersion parameter (a type of hyper-parameter)

• tidy() One-line summary

• set_n_draw() Specify number of prior or posterior draws

Forecast

• forecast() Use model to obtain estimates for future periods

Check model

• replicate_data() Compare real and simulated data

• report_sim() Simulation study of model

Data

• datasets Overview of datasets

• svds Overview of scaled SVDs

Author(s)

Maintainer: John Bryant [email protected]

Authors:

• John Bryant [email protected]

• Junni Zhang [email protected]

Other contributors:

• Bayesian Demography Limited [copyright holder]

See Also

Useful links:

• https://bayesiandemography.github.io/bage/

• https://github.com/bayesiandemography/bage

• Report bugs at https://github.com/bayesiandemography/bage/issues

---

Autoregressive Prior

Description

Use an autoregressive process to model a main effect, or use multiple autoregressive processes to model an interaction. Typically used with time effects or with interactions that involve time.

Usage

AR(
  n_coef = 2,
  s = 1,
  shape1 = 5,
  shape2 = 5,
  along = NULL,
  con = c("none", "by")
)

Arguments

n_coef	Number of lagged terms in the model, ie the order of the model. Default is 2.
s	Scale for the prior for the innovations. Default is 1.
shape1, shape2	Parameters for beta-distribution prior for coefficients. Defaults are 5 and 5.
along	Name of the variable to be used as the 'along' variable. Only used with interactions.
con	Constraints on parameters. Current choices are "none" and "by". Default is "none". See below for details.

Details

If AR() is used with an interaction, then separate AR processes are constructed along the 'along' variable, within each combination of the 'by' variables.

By default, the autoregressive processes have order 2. Alternative choices can be specified through the n_coef argument.

Argument s controls the size of innovations. Smaller values for s tend to give smoother estimates.

Value

An object of class "bage_prior_ar".

Mathematical details

When AR() is used with a main effect,

$\beta_j = \phi_1 \beta_{j-1} + \cdots + \phi_{\mathtt{n\_coef}} \beta_{j-\mathtt{n\_coef}} + \epsilon_j$

$\epsilon_j \sim \text{N}(0, \omega^2),$

and when it is used with an interaction,

$\beta_{u,v} = \phi_1 \beta_{u,v-1} + \cdots + \phi_{\mathtt{n\_coef}} \beta_{u,v-\mathtt{n\_coef}} + \epsilon_{u,v}$

$\epsilon_{u,v} \sim \text{N}(0, \omega^2),$

where

• $\pmb{\beta}$ is the main effect or interaction;

• $j$ denotes position within the main effect;

• $u$ denotes position within the 'by' variable(s) of the interaction; and

• $v$ denotes position within the 'along' variable of the interaction.

Internally, AR() derives a value for $\omega$ that gives every element of $\beta$ a marginal variance of $\tau^2$. Parameter $\tau$ has a half-normal prior

$\tau \sim \text{N}^+(0, \mathtt{s}^2).$

The correlation coefficients $\phi_1, \cdots, \phi_{\mathtt{n\_coef}}$ each have prior

$\phi_k \sim \text{Beta}(\mathtt{shape1}, \mathtt{shape2}).$

Constraints

With some combinations of terms and priors, the values of the intercept, main effects, and interactions are only weakly identified. This weak identifiability is typically harmless. However, in some applications, such as when trying to obtain interpretable values for main effects and interactions, it can be helpful to increase identifiability through the use of constraints, specified through the con argument.

Current options for con are:

• "none" No constraints. The default.

• "by" Only used in interaction terms that include 'along' and 'by' dimensions. Within each value of the 'along' dimension, terms across each 'by' dimension are constrained to sum to 0.

References

• AR() is based on the TMB function ARk

See Also

• AR1() Special case of AR(). Can be more numerically stable than higher-order models.

• Lin_AR(), Lin_AR1() Line with AR errors

• RW2_AR(), RW2_AR1() RW2 with AR errors

• priors Overview of priors implemented in bage

• set_prior() Specify prior for intercept, main effect, or interaction

• Mathematical Details vignette

Examples

AR(n_coef = 3)
AR(n_coef = 3, s = 2.4)
AR(along = "cohort")

---

Autoregressive Prior of Order 1

Description

Use an autoregressive process of order 1 to model a main effect, or use multiple AR1 processes to model an interaction. Typically used with time effects or with interactions that involve time.

Usage

AR1(
  s = 1,
  shape1 = 5,
  shape2 = 5,
  min = 0.8,
  max = 0.98,
  along = NULL,
  con = c("none", "by")
)

Arguments

s	Scale for the prior for the innovations. Default is 1.
shape1, shape2	Parameters for beta-distribution prior for coefficients. Defaults are 5 and 5.
min, max	Minimum and maximum values for autocorrelation coefficient. Defaults are 0.8 and 0.98.
along	Name of the variable to be used as the 'along' variable. Only used with interactions.
con	Constraints on parameters. Current choices are "none" and "by". Default is "none". See below for details.

Details

If AR() is used with an interaction, separate AR processes are constructed along the 'along' variable, within each combination of the 'by' variables.

Arguments min and max can be used to specify the permissible range for autocorrelation.

Argument s controls the size of innovations. Smaller values for s tend to give smoother estimates.

Value

An object of class "bage_prior_ar".

Mathematical details

When AR1() is used with a main effect,

$\beta_j = \phi \beta_{j-1} + \epsilon_j$

$\epsilon_j \sim \text{N}(0, \omega^2),$

and when it is used with an interaction,

$\beta_{u,v} = \phi \beta_{u,v-1} + \epsilon_{u,v}$

$\epsilon_{u,v} \sim \text{N}(0, \omega^2),$

where

• $\pmb{\beta}$ is the main effect or interaction;

• $j$ denotes position within the main effect;

• $u$ denotes position within the 'by' variable(s) of the interaction; and

• $v$ denotes position within the 'along' variable of the interaction.

Internally, AR1() derives a value for $\omega$ that gives every element of $\beta$ a marginal variance of $\tau^2$. Parameter $\tau$ has a half-normal prior

$\tau \sim \text{N}^+(0, \mathtt{s}^2),$

where s is provided by the user.

Coefficient $\phi$ is constrained to lie between min and max. Its prior distribution is

$\phi = (\mathtt{max} - \mathtt{min}) \phi' - \mathtt{min}$

where

$\phi' \sim \text{Beta}(\mathtt{shape1}, \mathtt{shape2}).$

Constraints

With some combinations of terms and priors, the values of the intercept, main effects, and interactions are only weakly identified. This weak identifiability is typically harmless. However, in some applications, such as when trying to obtain interpretable values for main effects and interactions, it can be helpful to increase identifiability through the use of constraints, specified through the con argument.

Current options for con are:

• "none" No constraints. The default.

• "by" Only used in interaction terms that include 'along' and 'by' dimensions. Within each value of the 'along' dimension, terms across each 'by' dimension are constrained to sum to 0.

References

• AR1() is based on the TMB function AR1

• The defaults for min and max are based on the defaults for forecast::ets().

See Also

• AR() Generalization of AR1()

• Lin_AR(), Lin_AR1() Line with AR errors

• RW2_AR(), RW2_AR1() RW2 with AR errors

• priors Overview of priors implemented in bage

• set_prior() Specify prior for intercept, main effect, or interaction

• Mathematical Details vignette

Examples

AR1()
AR1(min = 0, max = 1, s = 2.4)
AR1(along = "cohort")

---

Extract Data and Modeled Values

Description

Extract data and rates, probabilities, or means from a model object. The return value consists of the original data and one or more columns of modeled values.

Usage

## S3 method for class 'bage_mod'
augment(x, rows = NULL, quiet = FALSE, ...)

Arguments

x	Object of class "bage_mod", typically created with mod_pois(), mod_binom(), or mod_norm().
rows	Results are returned only for these rows of data. A logical vector, an index vector, or a logical expression evaluated within data. Optional
quiet	Whether to suppress messages. Default is FALSE.
...	Unused. Included for generic consistency only.

Details

The rows argument can be used to obtain results for a subset within data. This is faster, and uses less memory, than generating results for the whole dataset and then subsetting.

Value

A tibble, with the original data plus one or more of the following columns:

• ⁠.<outcome>⁠ Corrected or extended version of the outcome variable, in applications where the outcome variable has missing values, or a data model is being used.

• .observed 'Direct' estimates of rates or probabilities, ie counts divided by exposure or size (in Poisson and binomial models.)

• .fitted Draws of rates, probabilities, or means.

• .expected Draws of expected values for rates or probabilities (in Poisson that include exposure, or in binomial models.)

Uncertain quantities are represented using rvecs.

Fitted vs unfitted models

augment() is typically called on a fitted model. In this case, the modeled values are draws from the joint posterior distribution for rates, probabilities, or means.

augment() can, however, be called on an unfitted model. In this case, the modeled values are draws from the joint prior distribution. In other words, the modeled values are informed by model priors, and by values for exposure, size, or weights, but not by observed outcomes.

Imputed values for outcome variable

augment() automatically imputes any missing values for the outcome variable. If outcome variable var has one or more NAs, then augment creates a variable .var holding original and imputed values.

Data model for outcome variable

If the overall model includes a data model for the outcome variable var, then augment() creates a new variable .var containing estimates of the true value for the outcome.

See Also

• components() Extract values for hyper-parameters

• dispersion() Extract values for dispersion

• tidy() Short summary of a model

• mod_pois() Specify a Poisson model

• mod_binom() Specify a binomial model

• mod_norm() Specify a normal model

• fit() Fit a model

• is_fitted() See if a model has been fitted

• unfit() Reset a model

• datamods Overview of data models implemented in bage

Examples

set.seed(0)

## specify model
mod <- mod_pois(divorces ~ age + sex + time,
                data = nzl_divorces,
                exposure = population) |>
  set_n_draw(n_draw = 100) ## smaller sample, so 'augment' faster

## fit model
mod <- mod |>
  fit()

## draw from the posterior distribution
mod |>
  augment()

## results for females only
mod |>
  augment(rows = sex == "Female")

## insert a missing value into outcome variable
divorces_missing <- nzl_divorces
divorces_missing$divorces[1] <- NA

## fitting model and calling 'augument'A
## creates a new variable called '.divorces'
## holding observed and imputed values
mod_pois(divorces ~ age + sex + time,
         data = divorces_missing,
         exposure = population) |>
  fit() |>
  augment()

## specifying a data model for the
## original data also leads to a new
## variable called '.divorces'
mod_pois(divorces ~ age + sex + time,
         data = nzl_divorces,
         exposure = population) |>
  set_datamod_outcome_rr3() |>
  fit() |>
  augment()

---

Extract Values for Hyper-Parameters

Description

Extract values for hyper-parameters from a model object. Hyper-parameters include

• main effects and interactions,

• dispersion,

• trends, seasonal effects, errors,

• SVD, spline, and covariate coefficients,

• standard deviations, correlation coefficients.

Usage

## S3 method for class 'bage_mod'
components(object, quiet = FALSE, original_scale = FALSE, ...)

Arguments

object	Object of class "bage_mod", typically created with mod_pois(), mod_binom(), or mod_norm().
quiet	Whether to suppress messages. Default is FALSE.
original_scale	Whether values for "effect", "trend", "season", "error" and "disp" components from a normal model are on the original scale or the transformed scale. Default is FALSE.
...	Unused. Included for generic consistency only.

Value

A tibble with four columns columns:

The return value contains the following columns:

• term Model term that the hyper-parameter belongs to.

• component Component within term.

• level Element within component .

• .fitted An rvec containing draws from the posterior distribution.

Fitted vs unfitted models

components() is typically called on a fitted model. In this case, the values returned are draws from the joint posterior distribution for the hyper-parameters in the model.

components() can, however, be called on an unfitted model. In this case, the values returned are draws from the joint prior distribution. In other words, the values incorporate model priors, and any exposure, size, or weights argument, but not observed outcomes.

Scaling and Normal models

Internally, models created with mod_norm() are fitted using transformed versions of the outcome and weights variables. By default, when components() is used with these models, it returns values for .fitted that are based on the transformed versions. To instead obtain values for "effect", "trend", "season", "error" and "disp" that are based on the untransformed versions, set original_scale to TRUE.

See Also

• augment() Extract values for rates, means, or probabilities, together with original data

• dispersion() Extract values for dispersion

• tidy() Extract a one-line summary of a model

• mod_pois() Specify a Poisson model

• mod_binom() Specify a binomial model

• mod_norm() Specify a normal model

• fit() Fit a model

• is_fitted() See if a model has been fitted

• unfit() Reset a model

Examples

set.seed(0)

## specify model
mod <- mod_pois(injuries ~ age + sex + year,
                data = nzl_injuries,
                exposure = popn)

## extract prior distribution
## of hyper-parameters
mod |>
  components()

## fit model
mod <- mod |>
  fit()

## extract posterior distribution
## of hyper-parameters
mod |>
  components()

## fit normal model
mod <- mod_norm(value ~ age * diag + year,
                data = nld_expenditure,
                weights = 1) |>
  fit()

## dispersion (= standard deviation in normal model)
## on the transformed scale
mod |>
  components() |>
  subset(component == "disp")

## disperson on the original scale
mod |>
  components(original_scale = TRUE) |>
  subset(component == "disp")

---

Extract Components used by SVD Summary

Description

Extract the matrix and offset used by a scaled SVD summary of a demographic database.

Usage

## S3 method for class 'bage_ssvd'
components(
  object,
  v = NULL,
  n_comp = NULL,
  indep = NULL,
  age_labels = NULL,
  ...
)

Arguments

object	An object of class "bage_ssvd".
v	Version of scaled SVD components to use. If no value is suppled, the most recent version is used.
n_comp	The number of components. The default is half the total number of components of object.
indep	Whether to use independent or joint SVDs for each sex/gender, if the data contains a sex/gender variable. The default is to use independent SVDs. To obtain results for the total population when the data contains a sex/gender variable, set indep to NA.
age_labels	Age labels for the desired age or age-sex profile. If no labels are supplied, the most detailed profile available is used.
...	Not currently used.

Value

A tibble with the offset and components.

Scaled SVDs of demographic databases in bage

• HMD Mortality rates from the Human Mortality Database.

• HFD Fertility rates from the Human Fertility Database.

• LFP Labor forcce participation rates from the OECD.

See Also

• generate() Randomly generate age-profiles, or age-sex profiles, based on a scaled SVD summary.

• SVD() SVD prior for terms involving age.

• SVD_AR1(), SVD_AR(), SVD_RW(), SVD_RW2() Dynamic SVD priors for terms involving age and time.

• poputils::age_labels() Generate age labels.

Examples

## females and males modeled independently
components(LFP, n_comp = 3)

## joint model for females and males
components(LFP, indep = FALSE, n_comp = 3)

## females and males combined
components(LFP, indep = NA, n_comp = 3)

## specify age groups
labels <- poputils::age_labels(type = "five", min = 15, max = 60)
components(LFP, age_labels = labels)

---

Information on Computations Performed During Model Fitting

Description

Get information on computations performed by function fit(). The information includes the total time used for fitting, and the time used for two particular tasks that can be slow: running the optimizer stats::nlminb(), and drawing from the multivariate normal returned by the TMB. It also includes values returned by the optimizer: the number of iterations needed, and messages about convergence.

Usage

computations(object)

Arguments

object

A fitted object of class "bage_mod".

Value

A tibble with the following variables:

• time_total Seconds used for whole fitting process.

• time_max Seconds used for optimisiation.

• time_draw Seconds used by function TMB::sdreport().

• iter Number of iterations required for optimization.

• message Message about convergence returned by optimizer.

See Also

• mod_pois(),mod_binom(),mod_norm() Specify a model

• fit() Fit a model

• tidy() Summarise a model

• set_n_draw() Specify number of posterior draws

Examples

mod <- mod_pois(divorces ~ age + sex + time,
                data = nzl_divorces,
                exposure = population) |>
  fit()

computations(mod)

---

Confidentialization

Description

The models for rates, probabilities, or means created with functions mod_pois(), mod_binom(), and mod_norm() can be extended by adding descriptions of confidentalization procedures applied to the outcome variable.

Details

Data models for outcome variable

Function	Confidentialization procedure	pois	binom	norm
set_confidential_rr3()	Outcome randomly rounded to multiple of 3	Yes	Yes	No

---

Scaled SVD Components from Census School Attendance Data

Description

An object of class "bage_ssvd" holding scaled SVD components derived from census data on school attendance. The attendance data is assembed by the United Nations Statistics Division. CSA holds 5 components.

Usage

CSA

Format

Object of class "bage_ssvd".

Versions:

• "v2025" (default). Data downloaded on 2025-11-05

Warning

Compared other age-sex patterns for other demographic processes such as mortality, age-sex patterns for school attendance show substantial variation across populations. More components may be needed to obtain satisfactory models of age-sex patterns for school attendance than for other processes.

Source

Derived from data in the "Population 5 to 24 years of age by school attendance, sex and urban/rural residence" table from the Population Censuses' Datasets database assembled by the United Nations Statistics Division. Code to create CSA is in folder ‘data-raw/ssvd_csa’ in the source code for the bage package.

See Also

• Scaled SVDs Overview of scaled SVDs implemented in bage

• SVD() A prior based on a scaled SVD

---

Data to Create Scaled SVD Object Based on World Marriage Database

Description

A subset of the data needed to produce a scaled SVD object, derived from data from the World Marriage Database. The data is formatted using function data_ssvd_wmd() in package bssvd.

Usage

data_wmd

Format

A tibble with 6 rows and with columns version, type, labels_age, labels_sexgender, matrix, and offset.

Source

Derived from data from the World Marriage Data 2019 database available on the United Nations Population Division website, and assembled by the UNPD from national census and survey data.

See Also

• ssvd() Function to create scaled SVD objects

• WMD_C Scaled SVD object based on a full set of World Marriage Database data.

• Scaled SVDs Overview of scaled SVDs implemented in bage

---

Data Models

Description

The models for rates, probabilities, or means created with functions mod_pois(), mod_binom(), and mod_norm() can be extended by adding data models, also referred to as measurement error models.

Details

Function	Assumptions about measurement error	Poisson	Binomial	Normal
set_datamod_miscount()	Reported outcome has undercount and overcount	Yes	No	No
set_datamod_undercount()	Reported outcome has undercount	Yes	Yes	No
set_datamod_overcount()	Reported outcome has overcount	Yes	No	No
set_datamod_noise()	Reported outcome unbiased, but with positive and negative measurement errors	Yes*	No	Yes
set_datamod_exposure()	Reported exposure unbiased, but with positive and negative measurement errors	Yes*	No	No

*Models with no dispersion term for rates.

---

Datasets

Description

Datasets in bage

Details

dataset	Outcome	Variables	Country
isl_deaths	Deaths	age, sex, time, deaths, popn	Iceland
kor_births	Births	age, region, time, popn, gdp_pc_2023, dens2020	South Korea
nld_expenditure	Health expenditure	diag, age, year, value	Netherlands
nzl_divorces	Divorces	age, sex, time, divorces, population	New Zealand
nzl_households	One-person households	age, region, year, oneperson, total	New Zealand
nzl_injuries	Accidental deaths	age, sex, ethnicity, year, injuries, popn	New Zealand
prt_deaths	Deaths	age, time, deaths, exposure	Portugal
swe_infant	Infant deaths	county, time, births, deaths	Sweden
usa_deaths	Accidental deaths	month, deaths	United States

---

Extract Values for Dispersion

Description

Extract values for the 'dispersion' parameter from a model object.

Usage

dispersion(object, quiet = FALSE, original_scale = FALSE)

Arguments

object	Object of class "bage_mod", typically created with mod_pois(), mod_binom(), or mod_norm().
quiet	Whether to suppress messages. Default is FALSE.
original_scale	Whether values for disperson are on the original scale or the transformed scale. Default is FALSE.

Value

An rvec (or NULL if the model does not include a dispersion parameter.)

Fitted vs unfitted models

dispersion() is typically called on a fitted model. In this case, the values for dispersion are draws from the posterior distribution. dispersion() can, however, be called on an unfitted model. In this case, the values are drawn from the prior distribution.

Scaling and Normal models

Internally, models created with mod_norm() are fitted using transformed versions of the outcome and weights variables. By default, when dispersion() is used with these models, it returns values on the transformed scale. To instead obtain values on the untransformed scale, set original_scale to TRUE.

See Also

• components() Extract values for hyper-parameters, including dispersion

• set_disp() Specify a prior for dispersion

Examples

set.seed(0)

## specify model
mod <- mod_pois(injuries ~ age + sex + year,
                data = nzl_injuries,
                exposure = popn)

## prior distribution
mod |>
  dispersion()

## fit model
mod <- mod |>
  fit()

## posterior distribution
mod |>
  dispersion()

## fit normal model
mod <- mod_norm(value ~ age * diag + year,
                data = nld_expenditure,
                weights = 1) |>
  fit()

## values on the transformed scale
mod |>
  dispersion()

## values on the original scale
mod |>
  dispersion(original_scale = TRUE)

---

Damped Random Walk Prior

Description

Use a damped random walk as a model for a main effect, or use multiple damped random walks as a model for an interaction. Typically used with terms that involve time, particularly in forecasts. Damping often improves forecast accuracy.

Usage

DRW(
  s = 1,
  sd = 1,
  shape1 = 5,
  shape2 = 5,
  min = 0.8,
  max = 0.98,
  along = NULL,
  con = c("none", "by")
)

Arguments

s	Scale for the prior for the innovations. Default is 1.
sd	Standard deviation of initial value. Default is 1. Can be 0.
shape1, shape2	Parameters for beta-distribution prior for damping coefficient. Defaults are 5 and 5.
min, max	Minimum and maximum values for damping coefficient. Defaults are 0.8 and 0.98.
along	Name of the variable to be used as the 'along' variable. Only used with interactions.
con	Constraints on parameters. Current choices are "none" and "by". Default is "none". See below for details.

Details

If DRW() is used with an interaction, a separate damped random walk is constructed within each combination of the 'by' variables.

Arguments min and max can be used to control the amount of damping that occurs.

Argument s controls the size of innovations. Smaller values for s tend to produce smoother series.

Argument sd controls variance in initial values. Setting sd to 0 fixes initial values at 0.

Value

An object of class "bage_prior_drwrandom" or "bage_prior_drwzero".

Mathematical details

When DRW() is used with a main effect,

$\beta_1 \sim \text{N}(0, \mathtt{sd}^2)$

$\beta_j \sim \text{N}(\phi \beta_{j-1}, \tau^2), \quad j > 1$

and when it is used with an interaction,

$\beta_{u,1} \sim \text{N}(0, \mathtt{sd}^2)$

$\beta_{u,v} \sim \text{N}(\phi \beta_{u,v-1}, \tau^2), \quad v > 1$

where

• $\pmb{\beta}$ is the main effect or interaction;

• $\phi$ is the damping coefficient;

• $j$ denotes position within the main effect;

• $u$ denotes position within the 'by' variable(s) of the interaction; and

• $v$ denotes position within the 'along' variable of the interaction.

Coefficient $\phi$ is constrained to lie between min and max. Its prior distribution is

$\phi = (\mathtt{max} - \mathtt{min}) \phi' - \mathtt{min}$

where

$\phi' \sim \text{Beta}(\mathtt{shape1}, \mathtt{shape2}).$

Standard deviation $\tau$ has a half-normal prior

$\tau \sim \text{N}^+(0, \mathtt{s}^2),$

where s is provided by the user.

DRW() has the same basic structure as AR1(). However, in DRW(), $\tau$ controls the variance of the innovations, but in AR1() $\tau$ controls the marginal variance of each $\beta_j$ or $\beta_{u,v}$.

Constraints

With some combinations of terms and priors, the values of the intercept, main effects, and interactions are only weakly identified. This weak identifiability is typically harmless. However, in some applications, such as when trying to obtain interpretable values for main effects and interactions, it can be helpful to increase identifiability through the use of constraints, specified through the con argument.

Current options for con are:

• "none" No constraints. The default.

• "by" Only used in interaction terms that include 'along' and 'by' dimensions. Within each value of the 'along' dimension, terms across each 'by' dimension are constrained to sum to 0.

See Also

• DRW2() Damped second-order random walk

• RW() Random walk, without damping

• RW2() Second-order random walk, without damping

• RW_Seas() Random walk with seasonal effect

• AR() Autoregressive with order k

• AR1() Autoregressive with order 1

• Sp() Smoothing via splines

• SVD() Smoothing over age using singular value decomposition

• priors Overview of priors implemented in bage

• set_prior() Specify prior for intercept, main effect, or interaction

• Mathematical Details vignette

Examples

DRW()
DRW(min = 0, max = 1)
DRW(sd = 0)

---

Damped Second-Order Random Walk Prior

Description

Use a damped second-order random walk as a model for a main effect, or use multiple second-order random walks as a model for an interaction. A damped second-order random walk is a random walk with drift where the drift term varies, with a tendency to converge on zero. It is typically used with terms that involve time, where there are sustained trends upward or downward. Damping often improves forecast accuracy.

Usage

DRW2(
  s = 1,
  sd = 1,
  sd_slope = 1,
  shape1 = 5,
  shape2 = 5,
  min = 0.8,
  max = 0.98,
  along = NULL,
  con = c("none", "by")
)

Arguments

s	Scale for the prior for the innovations. Default is 1.
sd	Standard deviation of initial value. Default is 1. Can be 0.
sd_slope	Standard deviation of initial slope. Default is 1.
shape1, shape2	Parameters for beta-distribution prior for damping coefficient. Defaults are 5 and 5.
min, max	Minimum and maximum values for damping coefficient. Defaults are 0.8 and 0.98.
along	Name of the variable to be used as the 'along' variable. Only used with interactions.
con	Constraints on parameters. Current choices are "none" and "by". Default is "none". See below for details.

Details

If DRW2() is used with an interaction, a separate damped random walk is constructed within each combination of the 'by' variables.

Arguments min and max can be used to control the amount of damping that occurs.

Argument s controls the size of innovations. Smaller values for s tend to give smoother series.

Argument sd controls variance in initial values. Setting sd to 0 fixes initial values at 0.

Argument sd_slope controls variance in the initial slope.

Value

An object of class "bage_prior_drw2random" or "bage_prior_drw2zero".

Mathematical details

When DRW2() is used with a main effect,

$\beta_1 \sim \text{N}(0, \mathtt{sd}^2)$

$\beta_2 \sim \text{N}(\beta_1, \mathtt{sd\_slope}^2)$

$\beta_j \sim \text{N}(\beta_{j-1} + \phi (\beta_{j-1} - \beta_{j-2}), \tau^2), \quad j = 2, \cdots, J$

and when it is used with an interaction,

$\beta_{u,1} \sim \text{N}(0, \mathtt{sd}^2)$

$\beta_{u,2} \sim \text{N}(\beta_{u,1}, \mathtt{sd\_slope}^2)$

$\beta_{u,v} \sim \text{N}(\beta_{u,v-1} + \phi (\beta_{u,v-1} - \beta_{u,v-2}), \tau^2), \quad v = 3, \cdots, V$

where

• $\pmb{\beta}$ is the main effect or interaction;

• $\phi$ is the damping coefficient;

• $j$ denotes position within the main effect;

• $u$ denotes position within the 'by' variable(s) of the interaction; and

• $v$ denotes position within the 'along' variable of the interaction.

Coefficient $\phi$ is constrained to lie between min and max. Its prior distribution is

$\phi = (\mathtt{max} - \mathtt{min}) \phi' - \mathtt{min}$

where

$\phi' \sim \text{Beta}(\mathtt{shape1}, \mathtt{shape2}).$

Standard deviation $\tau$ has a half-normal prior

$\tau \sim \text{N}^+(0, \mathtt{s}^2),$

where s is provided by the user.

Constraints

With some combinations of terms and priors, the values of the intercept, main effects, and interactions are only weakly identified. This weak identifiability is typically harmless. However, in some applications, such as when trying to obtain interpretable values for main effects and interactions, it can be helpful to increase identifiability through the use of constraints, specified through the con argument.

Current options for con are:

• "none" No constraints. The default.

• "by" Only used in interaction terms that include 'along' and 'by' dimensions. Within each value of the 'along' dimension, terms across each 'by' dimension are constrained to sum to 0.

See Also

• DRW() Damped first-order random walk

• RW2() Second-order random walk, without damping

• RW2_Seas() Second order random walk with seasonal effect

• AR() Autoregressive with order k

• AR1() Autoregressive with order 1

• Sp() Smoothing via splines

• SVD() Smoothing over age via singular value decomposition

• priors Overview of priors implemented in bage

• set_prior() Specify prior for intercept, main effect, or interaction

• Mathematical Details vignette

Examples

DRW2()
DRW2(s = 0.5)
DRW2(min = 0, max = 1)

---

Fit a Model

Description

Derive the posterior distribution for a model.

Usage

## S3 method for class 'bage_mod'
fit(
  object,
  method = c("standard", "inner-outer"),
  vars_inner = NULL,
  optimizer = c("multi", "nlminb", "BFGS", "CG"),
  quiet = TRUE,
  max_jitter = 1e-04,
  start_oldpar = FALSE,
  ...
)

Arguments

object	A bage_mod object, created with mod_pois(), mod_binom(), or mod_norm().
method	Estimation method. Current choices are "standard" (the default) and "inner-outer". See below for details.
vars_inner	Names of variables to use for inner model when method is ⁠"inner-outer". If ⁠NULL⁠(the default)⁠vars_inner' is the age, sex/gender, and time variables.
optimizer	Which optimizer to use. Current choices are "multi", "nlminb", "BFGS", and "CG". Default is "multi". See below for details.
quiet	Whether to suppress messages from optimizer. Default is TRUE.
max_jitter	Maximum quantity to add to diagonal of precision matrix, if Cholesky factorization is failing. Default is 0.0001.
start_oldpar	Whether the optimizer should start at previous estimates. Used only when fit() is being called on a fitted model. Default is FALSE.
...	Not currently used.

Value

A bage_mod object

Estimation methods

When method is "standard" (the default), all parameters, other than the lowest-level rates, probabilities, or means are jointly estimated within TMB.

When method is "inner-outer", estimation is carried out in multiple steps, which, in large models, can sometimes reduce computation times. In Step 1, a model only using the inner variables is fitted to the data. In Step 2, a model only using the outer variables is fitted to the data. In Step 3, values for dispersion are calculated. Parameter estimates from steps 1, 2, and 3 are then combined.

Optimizer

The choices for the optimizer argument are:

• "multi" Try "nlminb", and if that fails, restart from the parameter values where "nlminb" stopped, using "BFGS". The default.

• "nlminb" stats::nlminb()

• "BFGS" stats::optim() using method "BFGS".

• "GC" stats::optim() using method "CG" (conjugate gradient).

Cholesky factorization and max_jitter

Sampling from the posterior distribution requires performing a Cholesky factorization of the precision matrix returned by TMB. This factorization sometimes fails because of numerical problems. Adding a small quantity to the diagonal of the precision matrix can alleviate numerical problems, while potentially reducing accuracy. If the Cholesky factorization initially fails, bage will try again with progressively larger quantities added to the diagonal, up to the maximum set by max_jitter. Increasing the value of max_jitter can help suppress numerical problems. A safer strategy, however, is to simplify the model, or to use more informative priors.

Aggregation

Up to version 0.9.8 of bage, fit() always aggregated across cells with identical values of the predictor variables in formula (ie the variables to the right of ~) before fitting. For instance, if a dataset contained deaths and population disaggregated by age and sex, but the model formula was deaths ~ age, then fit() would aggregate deaths and population within each age category before fitting the model. From version 0.9.9, fit() only aggregates across cells with identical values if no data model is used, and if the model is Poisson with dispersion set to 0 or is normal. Note that this change in behavior has no effect on most models, since most models include all variables used to classify outcomes.

See Also

• mod_pois() Specify a Poisson model

• mod_binom() Specify a binomial model

• mod_norm() Specify a normal model

• augment() Extract values for rates, probabilities, or means, together with original data

• components() Extract values for hyper-parameters

• dispersion() Extract values for dispersion

• forecast() Forecast, based on a model

• report_sim() Simulation study of a model

• unfit() Reset a model

• is_fitted() Check if a model has been fitted

• Mathematical Details vignette

Examples

## specify model
mod <- mod_pois(injuries ~ age + sex + year,
                data = nzl_injuries,
                exposure = popn)

## examine unfitted model
mod

## fit model
mod <- fit(mod)

## examine fitted model
mod

## extract rates
aug <- augment(mod)
aug

## extract hyper-parameters
comp <- components(mod)
comp

---

Use a Model to Make a Forecast

Description

Forecast rates, probabilities, means, and other model parameters.

Usage

## S3 method for class 'bage_mod'
forecast(
  object,
  newdata = NULL,
  labels = NULL,
  output = c("augment", "components"),
  include_estimates = FALSE,
  rows = NULL,
  quiet = FALSE,
  ...
)

Arguments

object	A bage_mod object, typically created with mod_pois(), mod_binom(), or mod_norm().
newdata	Data frame with data for future periods.
labels	Labels for future time periods.
output	Type of output returned. "augment" (the default) or "components".
include_estimates	Whether to include historical estimates along with the forecasts. Default is FALSE.
rows	A logical vector, an index vector, or a logical expression specifying which rows to return from an "augment" forecast. Can only be used if include_estimates is FALSE.
quiet	Whether to suppress messages. Default is FALSE.
...	Not currently used.

Value

A tibble.

How the forecasts are constructed

Internally, the steps involved in a forecast are:

1. Forecast time-varying main effects and interactions, e.g. a time main effect, or an age-time interaction.

2. Combine forecasts for the time-varying main effects and interactions with non-time-varying parameters, e.g. age effects or dispersion.

3. Use the combined parameters to generate values for rates, probabilities or means.

4. Optionally, generate values for the outcome variable.

forecast() generates values for the outcome variable when,

• output is "augment",

• a value has been supplied for newdata,

• newdata included a value for the exposure, size, or weights variable (except if exposure = 1 or weights = 1 in the original call to mod_pois() or mod_norm()).

Mathematical Details gives more details on the internal calculations in forecasting.

Output format

When output is "augment" (the default), the return value from forecast() looks like output from function augment(). When output is "components", the return value looks like output from components().

When include_estimates is FALSE (the default), the output of forecast() excludes values for time-varying parameters for the period covered by the data. When include_estimates is TRUE, the output includes these values. Setting include_estimates to TRUE can be helpful when creating graphs that combine estimates and forecasts.

Forecasting with covariates

Models that contain covariates can be used in forecasts, provided that

• all coefficients (the $\zeta_p$) are estimated from historical data via fit(), and

• if any covariates (the columns of $\pmb{Z}$) are time-varying, then future values for these covariates are supplied via the newdata argument.

Forecasting with data models

Models that contain data models can be used in forecasts, provided that

• the data models have no time-varying parameters, or

• future values for time-varying parameters are supplied when the data model is first specified.

For examples, see the Data Models vignette.

Fitted and unfitted models

forecast() is typically used with a fitted model, i.e. a model in which parameter values have been estimated from the data. The resulting forecasts reflect data and priors.

forecast() can, however, be used with an unfitted model. In this case, the forecasts are based entirely on the priors. See below for an example. Experimenting with forecasts based entirely on the priors can be helpful for choosing an appropriate model.

Warning

The interface for forecast() has not been finalised.

See Also

• mod_pois() Specify a Poisson model

• mod_binom() Specify a binomial model

• mod_norm() Specify a normal model

• fit() Fit a model

• augment() Extract values for rates, probabilities, or means, together with original data

• components() Extract values for hyper-parameters

• Mathematical Details vignette

Examples

## specify and fit model
mod <- mod_pois(injuries ~ age * sex + ethnicity + year,
                data = nzl_injuries,
                exposure = popn) |>
  fit()
mod

## forecasts
mod |>
  forecast(labels = 2019:2024)

## combined estimates and forecasts
mod |>
  forecast(labels = 2019:2024,
           include_estimates = TRUE)

## only selected rows
mod |>
  forecast(labels = 2019:2024,
           rows = age == "40-44")

## hyper-parameters
mod |>
  forecast(labels = 2019:2024,
           output = "components")

## hold back some data and forecast
library(dplyr, warn.conflicts = FALSE)
data_historical <- nzl_injuries |>
  filter(year <= 2015)
data_forecast <- nzl_injuries |>
  filter(year > 2015) |>
  mutate(injuries = NA)
mod_pois(injuries ~ age * sex + ethnicity + year,
         data = data_historical,
         exposure = popn) |>
  fit() |>
  forecast(newdata = data_forecast)

## forecast using GDP per capita in 2023 as a covariate
mod_births <- mod_pois(births ~ age * region + time,
                       data = kor_births,
                       exposure = popn) |>
  set_covariates(~ gdp_pc_2023) |>
  fit()
mod_births |>
  forecast(labels = 2024:2025)

---

Generate Values from Priors

Description

Generate draws from priors for model terms.

Usage

## S3 method for class 'bage_prior_ar'
generate(x, n_along = 20, n_by = 1, n_draw = 25, ...)

## S3 method for class 'bage_prior_drwrandom'
generate(x, n_along = 20, n_by = 1, n_draw = 25, ...)

## S3 method for class 'bage_prior_drwzero'
generate(x, n_along = 20, n_by = 1, n_draw = 25, ...)

## S3 method for class 'bage_prior_drw2random'
generate(x, n_along = 20, n_by = 1, n_draw = 25, ...)

## S3 method for class 'bage_prior_drw2zero'
generate(x, n_along = 20, n_by = 1, n_draw = 25, ...)

## S3 method for class 'bage_prior_known'
generate(x, n_element = 20, n_draw = 25, ...)

## S3 method for class 'bage_prior_lin'
generate(x, n_along = 20, n_by = 1, n_draw = 25, ...)

## S3 method for class 'bage_prior_linar'
generate(x, n_along = 20, n_by = 1, n_draw = 25, ...)

## S3 method for class 'bage_prior_linex'
generate(x, n_along = 20, n_by = 1, n_draw = 25, ...)

## S3 method for class 'bage_prior_norm'
generate(x, n_element = 20, n_draw = 25, ...)

## S3 method for class 'bage_prior_normfixed'
generate(x, n_element = 20, n_draw = 25, ...)

## S3 method for class 'bage_prior_rwrandom'
generate(x, n_along = 20, n_by = 1, n_draw = 25, ...)

## S3 method for class 'bage_prior_rwrandomseasfix'
generate(x, n_along = 20, n_by = 1, n_draw = 25, ...)

## S3 method for class 'bage_prior_rwrandomseasvary'
generate(x, n_along = 20, n_by = 1, n_draw = 25, ...)

## S3 method for class 'bage_prior_rwzero'
generate(x, n_along = 20, n_by = 1, n_draw = 25, ...)

## S3 method for class 'bage_prior_rwzeroseasfix'
generate(x, n_along = 20, n_by = 1, n_draw = 25, ...)

## S3 method for class 'bage_prior_rwzeroseasvary'
generate(x, n_along = 20, n_by = 1, n_draw = 25, ...)

## S3 method for class 'bage_prior_rw2random'
generate(x, n_along = 20, n_by = 1, n_draw = 25, ...)

## S3 method for class 'bage_prior_rw2randomar'
generate(x, n_along = 20, n_by = 1, n_draw = 25, ...)

## S3 method for class 'bage_prior_rw2randomseasfix'
generate(x, n_along = 20, n_by = 1, n_draw = 25, ...)

## S3 method for class 'bage_prior_rw2randomseasvary'
generate(x, n_along = 20, n_by = 1, n_draw = 25, ...)

## S3 method for class 'bage_prior_rw2zero'
generate(x, n_along = 20, n_by = 1, n_draw = 25, ...)

## S3 method for class 'bage_prior_rw2zeroar'
generate(x, n_along = 20, n_by = 1, n_draw = 25, ...)

## S3 method for class 'bage_prior_rw2zeroseasfix'
generate(x, n_along = 20, n_by = 1, n_draw = 25, ...)

## S3 method for class 'bage_prior_rw2zeroseasvary'
generate(x, n_along = 20, n_by = 1, n_draw = 25, ...)

## S3 method for class 'bage_prior_spline'
generate(x, n_along = 20, n_by = 1, n_draw = 25, ...)

## S3 method for class 'bage_prior_svd'
generate(x, n_element = 1, n_draw = 25, ...)

## S3 method for class 'bage_prior_svd_ar'
generate(x, n_along = 5, n_by = 1, n_draw = 25, ...)

## S3 method for class 'bage_prior_svd_drwrandom'
generate(x, n_along = 5, n_by = 1, n_draw = 25, ...)

## S3 method for class 'bage_prior_svd_drwzero'
generate(x, n_along = 5, n_by = 1, n_draw = 25, ...)

## S3 method for class 'bage_prior_svd_drw2random'
generate(x, n_along = 5, n_by = 1, n_draw = 25, ...)

## S3 method for class 'bage_prior_svd_drw2zero'
generate(x, n_along = 5, n_by = 1, n_draw = 25, ...)

## S3 method for class 'bage_prior_svd_lin'
generate(x, n_along = 5, n_by = 1, n_draw = 25, ...)

## S3 method for class 'bage_prior_svd_linex'
generate(x, n_along = 5, n_by = 1, n_draw = 25, ...)

## S3 method for class 'bage_prior_svd_rwrandom'
generate(x, n_along = 5, n_by = 1, n_draw = 25, ...)

## S3 method for class 'bage_prior_svd_rwzero'
generate(x, n_along = 5, n_by = 1, n_draw = 25, ...)

## S3 method for class 'bage_prior_svd_rw2random'
generate(x, n_along = 5, n_by = 1, n_draw = 25, ...)

## S3 method for class 'bage_prior_svd_rw2zero'
generate(x, n_along = 5, n_by = 1, n_draw = 25, ...)

Arguments

x	Object of class "bage_prior"
n_along	Number of elements of 'along' dimension. Default is 20.
n_by	Number of combinations of 'by' variables. Default is 1.
n_draw	Number of draws. Default is 25.
...	Unused. Included for generic consistency only.
n_element	Number of elements in term, in priors that do not distinguish 'along' and 'by' dimensions. Default is 20.

Details

Some priors distinguish between 'along' and 'by' dimensions, while others do not: see priors for a complete list. Arguments n_along and n_by are used with priors that make the distinction, and argument n_element is used with priors that do not.

Value

A tibble

See Also

• priors Overview of priors implemented in bage

Examples

## prior that distinguishes 'along' and 'by'
x <- RW()
generate(x, n_along = 10, n_by = 2)

## prior that does not distinguish
x <- N()
generate(x, n_element = 20)

## SVD_AR(), SVD_RW(), and SVD_RW2()
## distinguish 'along' and 'by'
x <- SVD_AR(HFD)
generate(x, n_along = 5, n_by = 2)

## SVD() does not
x <- SVD(HFD)
generate(x, n_element = 10)

---

Generate Random Age or Age-Sex Profiles

Description

Generate random age or age-sex profiles from an object of class "bage_ssvd". An object of class "bage_ssvd" holds results from an SVD decomposition of demographic data.

Usage

## S3 method for class 'bage_ssvd'
generate(
  x,
  v = NULL,
  n_draw = 20,
  n_comp = NULL,
  indep = NULL,
  age_labels = NULL,
  ...
)

Arguments

x	An object of class "bage_ssvd".
v	Version of data to use.
n_draw	Number of random draws to generate.
n_comp	The number of components. The default is half the total number of components of object.
indep	Whether to use independent or joint SVDs for each sex/gender, if the data contains a sex/gender variable. The default is to use independent SVDs. To obtain results for the total population when the data contains a sex/gender variable, set indep to NA.
age_labels	Age labels for the desired age or age-sex profile. If no labels are supplied, the most detailed profile available is used.
...	Unused. Included for generic consistency only.

Value

A tibble

Scaled SVDs of demographic databases in bage

• HMD Mortality rates from the Human Mortality Database.

• HFD Fertility rates from the Human Fertility Database.

• LFP Labor forcce participation rates from the OECD.

See Also

• components() Components used by SVD prior.

• SVD() SVD prior for term involving age.

• SVD_AR1(), SVD_AR(), SVD_RW(), SVD_RW2() Dynamic SVD priors for terms involving age and time.

• poputils::age_labels() Generate age labels.

Examples

## females and males modeled independently
generate(HMD) 

## joint model for females and males
generate(HMD, indep = FALSE) 

## SVD for females and males combined
generate(HMD, indep = NA)

## specify age groups
labels <- poputils::age_labels(type = "lt", max = 60)
generate(HMD, age_labels = labels)

---

Scaled SVD Components from Human Fertility Database

Description

An object of class "bage_ssvd" holding scaled SVD components derived from data from the Human Fertility Database. HFD holds 5 components.

Usage

HFD

Format

Object of class "bage_ssvd".

Versions:

• "v2025" (default) Data published on 2025-07-24

• "v2024" Data published on October 2024-10-23

Source

Derived from data from the Human Fertility Database. Max Planck Institute for Demographic Research (Germany) and Vienna Institute of Demography (Austria). Code to create HFD is in folder ‘data-raw/ssvd_hfd’ in the source code for the bage package.

See Also

• Scaled SVDs Overview of scaled SVDs implemented in bage

• SVD() A prior based on a scaled SVD

---

Scaled SVD Components from Human Internal Migration Database

Description

Objects of class "bage_ssvd" holding scaled SVD components derived from data from the Human Internal Migration Database. HIMD_P1, HIMD_P5, and HIMD_R each hold 5 components

Usage

HIMD_R

HIMD_P1

HIMD_P5

Format

Object of class "bage_ssvd".

Versions:

• "v2024" (default) Data published on 2024-10-23

Details

• HIMD_P1 is derived from data on 1-year migration probabilities, ie the probability that a person will migrate during a time interval of 1 year.

• HIMD_P5 is derived from data on 5-year migration probabilities, ie the probability that a person will migrate during a time interval of 5 years.

• HIMD_R is derived from data on 1-year migration probabilities, using the formula $r = -\log(1 - p)$.

Source

Dyrting, S. (2024, October 23). Data from: Estimating Complete Migration Probabilities from Grouped Data. Retrieved from osf.io/vmrfk on 1 September 2025. Code to create HIMD_R, HIMD_P1 and HIMD_P5 is in folder ‘data-raw/ssvd_himd’ in the source code for the bage package.

See Also

• Scaled SVDs Overview of scaled SVDs implemented in bage

• SVD() A prior based on a scaled SVD

---

Scaled SVD Components from Human Mortality Database

Description

An object of class "bage_ssvd" holding scaled SVD components derived from data from the Human Mortality Database. HMD holds 5 components.

Usage

HMD

Format

Object of class "bage_ssvd".

Versions:

• "v2025" (default) Data published on 2025-09-25, all years

• "v2025-50" Data published on 2025-09-25, 1950 and later

• "v2024" Data published on 2024-02-26, all years

Source

Derived from data from the Human Mortality Database. Max Planck Institute for Demographic Research (Germany), University of California, Berkeley (USA), and French Institute for Demographic Studies (France). Code to create HMD is in folder ‘data-raw/ssvd_hmd’ in the source code for the bage package.

See Also

• Scaled SVDs Overview of scaled SVDs implemented in bage

• SVD() A prior based on a scaled SVD

---

Test Whether a Model has Been Fitted

Description

Test whether fit() has been called on a model object.

Usage

is_fitted(x)

Arguments

x	An object of class "bage_mod".

Value

TRUE or FALSE

See Also

• mod_pois(), mod_binom(), mod_norm() to specify a model

• fit() to fit a model

Examples

mod <- mod_pois(injuries ~ age + sex + year,
                data = nzl_injuries,
                exposure = popn)
is_fitted(mod)
mod <- fit(mod)
is_fitted(mod)

---

Deaths in Iceland

Description

Deaths and mid-year populations in Iceland, by age, sex, and calendar year.

Usage

isl_deaths

Format

A tibble with 5,300 rows and the following columns:

• age Single year of age, up to "105+"

• sex "Female" and "Male"

• time Calendar year, 1998-2022

• deaths Counts of deaths

• popn Mid-year population

Source

Tables "Deaths by municipalities, sex and age 1981-2022", and "Average annual population by municipality, age and sex 1998-2022 - Current municipalities", on the Statistics Iceland website. Data downloaded on 12 July 2023.

See Also

• datasets Overview of datasets in bage

---

Known Prior

Description

Treat an intercept, a main effect, or an interaction as fixed and known.

Usage

Known(values)

Arguments

values

A numeric vector

Value

An object of class "bage_prior_known".

See Also

• NFix() Prior where level unknown, but variability known.

• priors Overview of priors implemented in bage

• set_prior() Specify prior for intercept, main effect, or interaction

• Mathematical Details vignette

Examples

Known(-2.3)
Known(c(0.1, 2, -0.11))

---

Births in South Korea

Description

Births and mid-year population by age of mother, region, and calendar year, 2011-2023, plus regional data on GDP per capita and population density.

Usage

kor_births

Format

A tibble with 1,872 rows and the following columns:

• age Five-year age groups from ⁠"10-14" to ⁠"50-54"'

• region Administrative region

• time Calendar year, 2011-2023

• births Counts of births

• popn Mid-year population

• gdp_pc_2023 Regional GDP per capita in 2023

• dens_2020 Regional population density (people per km-squared) in 2020

Source

Tables "Live Births by Age Group of Mother, Sex and Birth Order for Provinces", and "Resident Population in Five-Year Age Groups", on the Korean Statistical Information Service website. Data downloaded on 24 September 2024. Data on GDP per capita and population density from Wikipedia https://w.wiki/DMFA, data downloaded on 8 March 2025, and https://w.wiki/DMF9, data downloaded on 8 March 2025.

See Also

• datasets Overview of datasets in bage

---

Scaled SVD Components from OECD Labor Force Participation Data

Description

An object of class "bage_ssvd" holding scaled SVD components derived from labor force participation data assembled by the OECD. LFP holds 5 components.

Usage

LFP

Format

Object of class "bage_ssvd".

Versions:

• "v2025" Data downloaded on 2025-10-17

Source

Derived from data in the "Labor Force Indicators" table of the OECD Data Explorer. Code to create LFS is in folder ‘data-raw/ssvd_lfp’ in the source code for the bage package.

See Also

• Scaled SVDs Overview of scaled SVDs implemented in bage

• SVD() A prior based on a scaled SVD

---

Linear Prior with Independent Normal Errors

Description

Use a line or lines with independent normal errors to model a main effect or interaction. Typically used with time.

Usage

Lin(s = 1, mean_slope = 0, sd_slope = 1, along = NULL, con = c("none", "by"))

Arguments

s	Scale for the prior for the errors. Default is 1. Can be 0.
mean_slope	Mean in prior for slope of line. Default is 0.
sd_slope	Standard deviation in prior for slope of line. Default is 1.
along	Name of the variable to be used as the 'along' variable. Only used with interactions.
con	Constraints on parameters. Current choices are "none" and "by". Default is "none". See below for details.

Details

If Lin() is used with an interaction, then separate lines are constructed along the 'along' variable, within each combination of the 'by' variables.

Argument s controls the size of the errors. Smaller values give smoother estimates. s can be zero, in which case errors are zero, and all values lie exactly on straight lines. This is clearly a simplification, but it allows the prior to be used with very large interactions.

Argument sd_slope controls the size of the slopes of the lines. Larger values can give more steeply sloped lines.

Value

An object of class "bage_prior_lin".

Mathematical details

When Lin() is used with a main effect,

$\beta_j = (j - (J+1)/2) \eta + \epsilon_j$

$\eta \sim \text{N}(\mathtt{mean\_slope}, \mathtt{sd\_slope}^2)$

$\epsilon_j \sim \text{N}(0, \tau^2),$

and when it is used with an interaction,

$\beta_{u,v} = (v - (V + 1)/2) \eta_u + \epsilon_{u,v}$

$\eta_u \sim \text{N}(\mathtt{mean\_slope}, \mathtt{sd\_slope}^2)$

$\epsilon_{u,v} \sim \text{N}(0, \tau^2),$

where

• $\pmb{\beta}$ is the main effect or interaction;

• $j$ denotes position within the main effect;

• $v$ denotes position within the 'along' variable of the interaction; and

• $u$ denotes position within the 'by' variable(s) of the interaction.

Parameter $\tau$ has a half-normal prior

$\tau \sim \text{N}^+(0, \mathtt{s}^2).$

When $\mathtt{s} = 0$, the model reduces to

$\beta_j = (j - (J+1)/2) \eta$

$\eta \sim \text{N}(\mathtt{mean\_slope}, \mathtt{sd\_slope}^2)$

or

$\beta_{u,v} = (v = (V + 1)/2) \eta_u$

$\eta_u \sim \text{N}(\mathtt{mean\_slope}, \mathtt{sd\_slope}^2)$

.

Constraints

With some combinations of terms and priors, the values of the intercept, main effects, and interactions are only weakly identified. This weak identifiability is typically harmless. However, in some applications, such as when trying to obtain interpretable values for main effects and interactions, it can be helpful to increase identifiability through the use of constraints, specified through the con argument.

Current options for con are:

• "none" No constraints. The default.

• "by" Only used in interaction terms that include 'along' and 'by' dimensions. Within each value of the 'along' dimension, terms across each 'by' dimension are constrained to sum to 0.

See Also

• Lin_AR() Linear with AR errors

• Lin_AR1() Linear with AR1 errors

• RW2() Second-order random walk

• priors Overview of priors implemented in bage

• set_prior() Specify prior for intercept, main effect, or interaction

• Mathematical Details vignette

Examples

Lin()
Lin(s = 0.5, sd_slope = 2)
Lin(s = 0)
Lin(along = "cohort")

---

Linear Prior with Autoregressive Errors

Description

Use a line or lines with autoregressive errors to model a main effect or interaction. Typically used with time.

Usage

Lin_AR(
  n_coef = 2,
  s = 1,
  shape1 = 5,
  shape2 = 5,
  mean_slope = 0,
  sd_slope = 1,
  along = NULL,
  con = c("none", "by")
)

Arguments

n_coef	Number of lagged terms in the model, ie the order of the model. Default is 2.
s	Scale for the innovations in the AR process. Default is 1.
shape1, shape2	Parameters for beta-distribution prior for coefficients in the AR process. Defaults are 5 and 5.
mean_slope	Mean in prior for slope of line. Default is 0.
sd_slope	Standard deviation in the prior for the slope of the line. Larger values imply steeper slopes. Default is 1.
along	Name of the variable to be used as the 'along' variable. Only used with interactions.
con	Constraints on parameters. Current choices are "none" and "by". Default is "none". See below for details.

Details

If Lin_AR() is used with an interaction, separate lines are constructed along the 'along' variable, within each combination of the 'by' variables.

The order of the autoregressive errors is controlled by the n_coef argument. The default is 2.

Argument s controls the size of the innovations. Smaller values tend to give smoother estimates.

Argument sd_slope controls the slopes of the lines. Larger values can give more steeply sloped lines.

Value

An object of class "bage_prior_linar".

Mathematical details

When Lin_AR() is used with a main effect,

$\beta_j = \alpha_j + \epsilon_j$

$\alpha_j = (j - (J+1)/2) \eta$

$\epsilon_j = \phi_1 \epsilon_{j-1} + \cdots + \phi_{\mathtt{n\_coef}} \epsilon_{j-\mathtt{n\_coef}} + \varepsilon_j$

$\varepsilon_j \sim \text{N}(0, \omega^2),$

and when it is used with an interaction,

$\beta_{u,v} = \alpha_{u,v} + \epsilon_{u,v}$

$\alpha_{u,v} = (v - (V + 1)/2) \eta_u$

$\epsilon_{u,v} = \phi_1 \epsilon_{u,v-1} + \cdots + \phi_{\mathtt{n\_coef}} \epsilon_{u,v-\mathtt{n\_coef}} + \varepsilon_{u,v},$

$\varepsilon_{u,v} \sim \text{N}(0, \omega^2).$

where

• $\pmb{\beta}$ is the main effect or interaction;

• $j$ denotes position within the main effect;

• $u$ denotes position within the 'by' variable(s) of the interaction; and

• $v$ denotes position within the 'along' variable of the interaction.

The slopes have priors

$\eta \sim \text{N}(\mathtt{mean\_slope}, \mathtt{sd\_slope}^2)$

and

$\eta_u \sim \text{N}(\mathtt{mean\_slope}, \mathtt{sd\_slope}^2).$

Internally, Lin_AR() derives a value for $\omega$ that gives $\epsilon_j$ or $\epsilon_{u,v}$ a marginal variance of $\tau^2$. Parameter $\tau$ has a half-normal prior

$\tau \sim \text{N}^+(0, \mathtt{s}^2).$

The correlation coefficients $\phi_1, \cdots, \phi_{\mathtt{n\_coef}}$ each have prior

$0.5 \phi_k - 0.5 \sim \text{Beta}(\mathtt{shape1}, \mathtt{shape2}).$

Constraints

With some combinations of terms and priors, the values of the intercept, main effects, and interactions are only weakly identified. This weak identifiability is typically harmless. However, in some applications, such as when trying to obtain interpretable values for main effects and interactions, it can be helpful to increase identifiability through the use of constraints, specified through the con argument.

Current options for con are:

• "none" No constraints. The default.

• "by" Only used in interaction terms that include 'along' and 'by' dimensions. Within each value of the 'along' dimension, terms across each 'by' dimension are constrained to sum to 0.

See Also

• Lin_AR1() Special case of Lin_AR()

• Lin() Line with independent normal errors

• AR() AR process with no line

• priors Overview of priors implemented in bage

• set_prior() Specify prior for intercept, main effect, or interaction

• Mathematical Details vignette

Examples

Lin_AR()
Lin_AR(n_coef = 3, s = 0.5, sd_slope = 2)

---

Linear Prior with Autoregressive Errors of Order 1

Description

Use a line or lines with AR1 errors to model a main effect or interaction. Typically used with time.

Usage

Lin_AR1(
  s = 1,
  shape1 = 5,
  shape2 = 5,
  min = 0.8,
  max = 0.98,
  mean_slope = 0,
  sd_slope = 1,
  along = NULL,
  con = c("none", "by")
)

Arguments

s	Scale for the innovations in the AR process. Default is 1.
shape1, shape2	Parameters for beta-distribution prior for coefficients in the AR process. Defaults are 5 and 5.
min, max	Minimum and maximum values for autocorrelation coefficient in the AR process. Defaults are 0.8 and 0.98.
mean_slope	Mean in prior for slope of line. Default is 0.
sd_slope	Standard deviation in the prior for the slope of the line. Larger values imply steeper slopes. Default is 1.
along	Name of the variable to be used as the 'along' variable. Only used with interactions.
con	Constraints on parameters. Current choices are "none" and "by". Default is "none". See below for details.

Details

If Lin_AR1() is used with an interaction, separate lines are constructed along the 'along' variable, within each combination of the 'by' variables.

Arguments min and max can be used to specify the permissible range for autocorrelation.

Argument s controls the size of the innovations. Smaller values tend to give smoother estimates.

Argument sd_slope controls the slopes of the lines. Larger values can give more steeply sloped lines.

Value

An object of class "bage_prior_linar".

Mathematical details

When Lin_AR1() is being used with a main effect,

$\beta_j = \alpha_j + \epsilon_j$

$\alpha_j = (j - (J+1)/2) \eta$

$\epsilon_j = \phi \epsilon_{j-1} + \varepsilon_j$

$\varepsilon \sim \text{N}(0, \omega^2),$

and when it is used with an interaction,

$\beta_{u,v} = \alpha_{u,v} + \epsilon_{u,v}$

$\alpha_{u,v} = (v - (V + 1)/2) \eta_u$

$\epsilon_{u,v} = \phi \epsilon_{u,v-1} + \varepsilon_{u,v}$

$\varepsilon_{u,v} \sim \text{N}(0, \omega^2),$

where

• $\pmb{\beta}$ is the main effect or interaction;

• $j$ denotes position within the main effect;

• $u$ denotes position within the 'by' variable(s) of the interaction; and

• $v$ denotes position within the 'along' variable of the interaction.

The slopes have priors

$\eta \sim \text{N}(\mathtt{mean\_slope}, \mathtt{sd\_slope}^2)$

and

$\eta_u \sim \text{N}(\mathtt{mean\_slope}, \mathtt{sd\_slope}^2).$

Internally, Lin_AR1() derives a value for $\omega$ that gives $\epsilon_j$ or $\epsilon_{u,v}$ a marginal variance of $\tau^2$. Parameter $\tau$ has a half-normal prior

$\tau \sim \text{N}^+(0, \mathtt{s}^2),$

where a value for s is provided by the user.

Coefficient $\phi$ is constrained to lie between min and max. Its prior distribution is

$\phi = (\mathtt{max} - \mathtt{min}) \phi' - \mathtt{min}$

where

$\phi' \sim \text{Beta}(\mathtt{shape1}, \mathtt{shape2}).$

Constraints

With some combinations of terms and priors, the values of the intercept, main effects, and interactions are only weakly identified. This weak identifiability is typically harmless. However, in some applications, such as when trying to obtain interpretable values for main effects and interactions, it can be helpful to increase identifiability through the use of constraints, specified through the con argument.

Current options for con are:

• "none" No constraints. The default.

• "by" Only used in interaction terms that include 'along' and 'by' dimensions. Within each value of the 'along' dimension, terms across each 'by' dimension are constrained to sum to 0.

References

• The defaults for min and max are based on the defaults for forecast::ets().

See Also

• Lin_AR() Generalization of Lin_AR1()

• Lin() Line with independent normal errors

• AR1() AR1 process with no line

• priors Overview of priors implemented in bage

• set_prior() Specify prior for intercept, main effect, or interaction

• Mathematical Details vignette

Examples

Lin_AR1()
Lin_AR1(min = 0, s = 0.5, sd_slope = 2)

---

Specify a Binomial Model

Description

Specify a model where the outcome is drawn from a binomial distribution.

Usage

mod_binom(formula, data, size)

Arguments

formula	An R formula, specifying the outcome and predictors.
data	A data frame containing the outcome and predictor variables, and the number of trials.
size	Name of the variable giving the number of trials (with or without quote marks.)

Details

The model is hierarchical. The probabilities in the binomial distribution are described by a prior model formed from dimensions such as age, sex, and time. The terms for these dimension themselves have models, as described in priors. These priors all have defaults, which depend on the type of term (eg an intercept, an age main effect, or an age-time interaction.)

Value

An object of class bage_mod.

Mathematical details

The likelihood is

$y_i \sim \text{binomial}(\gamma_i; w_i)$

where

• subscript $i$ identifies some combination of the the classifying variables, such as age, sex, and time;

• $y_i$ is a count, such of number of births, such as age, sex, and region;

• $\gamma_i$ is a probability of 'success'; and

• $w_i$ is the number of trials.

The probabilities $\gamma_i$ are assumed to be drawn a beta distribution

$y_i \sim \text{Beta}(\xi^{-1} \mu_i, \xi^{-1} (1 - \mu_i))$

where

• $\mu_i$ is the expected value for $\gamma_i$; and

• $\xi$ governs dispersion (ie variance.)

Expected value $\mu_i$ equals, on a logit scale, the sum of terms formed from classifying variables,

$\text{logit} \mu_i = \sum_{m=0}^{M} \beta_{j_i^m}^{(m)}$

where

• $\beta^{0}$ is an intercept;

• $\beta^{(m)}$, $m = 1, \dots, M$, is a main effect or interaction; and

• $j_i^m$ is the element of $\beta^{(m)}$ associated with cell $i$.

The $\beta^{(m)}$ are given priors, as described in priors.

$\xi$ has an exponential prior with mean 1. Non-default values for the mean can be specified with set_disp().

The model for $\mu_i$ can also include covariates, as described in set_covariates().

See Also

• mod_pois() Specify Poisson model

• mod_norm() Specify normal model

• set_prior() Specify non-default prior for term

• set_disp() Specify non-default prior for dispersion

• fit() Fit a model

• augment() Extract values for probabilities, together with original data

• components() Extract values for hyper-parameters

• forecast() Forecast parameters and outcomes

• report_sim() Check model using simulation study

• replicate_data() Check model using replicate data

• Mathematical Details Detailed descriptions of models

Examples

mod <- mod_binom(oneperson ~ age:region + age:year,
                 data = nzl_households,
                 size = total)

---

Specify a Normal Model

Description

Specify a model where the outcome is drawn from a normal distribution.

Usage

mod_norm(formula, data, weights = NULL)

Arguments

formula	An R formula, specifying the outcome and predictors.
data	A data frame containing outcome, predictor, and, optionally, weights variables.
weights	Name of the weights variable, a 1, or a formula. See below for details.

Details

The model is hierarchical. The means in the normal distribution are described by a prior model formed from dimensions such as age, sex, and time. The terms for these dimension themselves have models, as described in priors. These priors all have defaults, which depend on the type of term (eg an intercept, an age main effect, or an age-time interaction.)

Value

An object of class bage_mod_norm.

Scaling of outcome and weights

Internally, mod_norm() scales the outcome variable to have mean 0 and standard deviation 1, and scales the weights to have mean 1. This scaling allows mod_norm() to use the same menu of priors as mod_pois() and mod_binom().

augment() always returns values on the original scale, rather than the transformed scale.

components() by default returns values on the transformed scale. But if original_scale is TRUE, it returns some types of values on the original scale. See components() for details.

Specifying weights

There are three options for creating an unweighted model:

• do not supply a value for the weights variable;

• set weights equal to NULL; or

• set weights equal to 1, though this option is deprecated, and will eventually be removed.

To create a weighted model, supply the name of the weighting variable in data, quoted or unquoted.

Mathematical details

The likelihood is

$y_i \sim \text{N}(\gamma_i, w_i^{-1} \sigma^2)$

where

• subscript $i$ identifies some combination of the classifying variables, such as age, sex, and time,

• $y_i$ is the value of the outcome variable,

• $w_i$ is a weight.

In some applications, $w_i$ is set to 1 for all $i$.

Internally, bage works with standardized versions of $\gamma_i$ and $\sigma^2$:

$\mu_i = (\gamma_i - \bar{y}) / s$

$\xi^2 = \sigma^2 / (\bar{w} s^2)$

where

$\bar{y} = \sum_{i=1}^n y_i / n$

$s = \sqrt{\sum_{i=1}^n (y_i - \bar{y})^2 / (n-1)}$

$\bar{w} = \sum_{i=1}^n w_i / n$

Mean parameter $\mu_i$ is modelled as the sum of terms formed from classifying variables and covariates,

$\mu_i = \sum_{m=0}^{M} \beta_{j_i^m}^{(m)}$

where

• $\beta^{0}$ is an intercept;

• $\beta^{(m)}$, $m = 1, \dots, M$, is a main effect or interaction; and

• $j_i^m$ is the element of $\beta^{(m)}$ associated with cell $i$,

The $\beta^{(m)}$ are given priors, as described in priors.

$\xi$ has an exponential prior with mean 1. Non-default values for the mean can be specified with set_disp().

The model for $\mu_i$ can also include covariates, as described in set_covariates().

See Also

• mod_pois() Specify Poisson model

• mod_binom() Specify binomial model

• set_prior() Specify non-default prior for term

• set_disp() Specify non-default prior for standard deviation

• fit() Fit a model

• augment() Extract values for means, together with original data

• components() Extract values for hyper-parameters

• forecast() Forecast parameters and outcomes

• report_sim() Check model using a simulation study

• replicate_data() Check model using replicate data data for a model

• Mathematical Details Detailed description of models

Examples

## model without weights
mod <- mod_norm(value ~ diag:age + year,
                data = nld_expenditure)

## model with weights
nld_expenditure$wt <- sqrt(nld_expenditure$value)
mod <- mod_norm(value ~ diag:age + year,
                data = nld_expenditure,
                weights = wt)

---