5.1.1 ‘Along’ and ‘by’ dimensions

Each \(\pmb{\beta}^{(m)}\), \(m > 0\), can be a main effect, involving a single dimension, or an interaction, involving two dimensions. Some priors, when applied to an interaction, treat one dimension, referred to as the ‘along’ dimension, differently from the remaining dimensions, referred to as ‘by’ dimensions. A random walk prior (Section 5.4), for instance, consists of an independent random walk along the ‘along’ dimension, within each combination of the ‘by’ dimensions.

We use \(v = 1, \cdots, V_m\) to denote position within the ‘along’ dimension, and \(u = 1, \cdots, V_m\) to denote position within a classification formed by the ‘by’ dimensions. When there are no sum-to-zero constraints (see below), \(U_m = \prod_k d_k\) where \(d_k\) is the number of elements in the \(k\)th ‘by’ variable. When there are sum-to-zero constraints, \(U_m = \prod_k (d_k - 1)\).

If a prior involves an ‘along’ dimension but the user does not specify one, the procedure for choosing a dimension is as follows:

• if the term involves time, use the time dimension;
• otherwise, if the term involves age, use the age dimension;
• otherwise, raise an error asking the user to explicitly specify a dimension.

5.1.2 Sum-to-zero constraints

Priors that involve ‘along’ dimensions can optionally include sum-to-zero constraints. If these constrains are applied, then within each element \(v\) of the ‘along’ dimension, the sum of the \(\beta_j^{(m)}\) across each ‘by’ dimension is zero. For instance, if \(\pmb{\beta}^{(m)}\) is an interaction between time, region, and sex, with time as the ‘along’ variable, then within each combination of time and region, the values for females and males sum to zero, and within each combination of time and sex, the values for regions sum to zero.

Except in the case of dynamic SVD-based priors (eg Sections 5.15), the sum-to-zero constraints are implemented internally by drawing values within an unrestricted lower-dimensional space, and then transforming to the restricted higher-dimensional space. For instance, a random walk prior for a time-region interaction with \(R\) regions consists of \(R-1\) unrestricted random walks along time, which are converted into \(R\) random walks that sum to zero across region. Matrices for transforming between the unrestricted and restricted spaces are constructed using the QR decomposition, as described in Section 1.8.1 of Wood (2017). With dynamic SVD-based priors, we draw values for the SVD coefficients with no constraints, convert these to unconstrained values for \(\pmb{\beta}^{(m)}\), and then subtract means.

5.1.3 Algorithm for assigning default priors

• If \(\pmb{\beta}^{(m)}\) has one or two elements, assign \(\pmb{\beta}^{(m)}\) a fixed-normal prior (Section 5.3);
• otherwise, if \(\pmb{\beta}^{(m)}\) involves time, assign \(\pmb{\beta}^{(m)}\) a random walk prior (Section 5.4) along the time dimension;
• otherwise, if \(\pmb{\beta}^{(m)}\) involves age, assign \(\pmb{\beta}^{(m)}\) a random walk prior (Section 5.4) along the age dimension;
• otherwise, assign \(\pmb{\beta}^{(m)}\) a normal prior (Section 5.2)

The intercept term \(\pmb{\beta}^{(0)}\) can only be given a fixed-normal prior (Section 5.3) or a Known prior (Section 5.19).

5.2.1 Model

Exchangeable normal

\[\begin{align} \beta_j^{(m)} & \sim \text{N}\left(0, \tau_m^2 \right) \\ \tau_m & \sim \text{N}^+\left(0, A_{\tau}^{(m)2}\right) \end{align}\]

5.2.2 Contribution to posterior density

\[\begin{equation} \text{N}(\tau_m \mid 0, A_{\tau}^{(m)2}) \prod_{j=1}^{J_m} \text{N}(\beta_j^{(m)} \mid 0, \tau_m^2) \end{equation}\]

5.2.3 Forecasting

\[\begin{equation} \beta_{J_m+h+1}^{(m)} \sim \text{N}(0, \tau_m^2) \end{equation}\]

5.2.4 Code

N(s = 1)

• s is \(A_{\tau}^{(m)}\). Defaults to 1.

5.3.1 Model

Exchangeable normal, with fixed standard deviation

\[\begin{equation} \beta_j^{(m)} \sim \text{N}\left(0, A_{\beta}^{(m)2}\right) \end{equation}\]

5.3.2 Contribution to posterior density

\[\begin{equation} \prod_{j=1}^{J_m} \text{N}(\beta_j^{(m)} \mid 0, A_{\beta}^{(m)2}) \end{equation}\]

5.3.3 Forecasting

\[\begin{equation} \beta_{J_m+h+1}^{(m)} \sim \text{N}(0, A_{\beta}^{(m)2}) \end{equation}\]

5.3.4 Code

NFix(sd = 1)

• sd is \(A_{\tau}^{(m)}\). Defaults to 1.

5.4.1 Model

Random walk

\[\begin{align} \beta_{u,1}^{(m)} & \sim \text{N}\left(0, (A_0^{(m)})^2\right) \\ \beta_{u,v}^{(m)} & \sim \text{N}(\beta_{u,v-1}^{(m)}, \tau_m^2), \quad v = 2, \cdots, V_m \\ \tau_m & \sim \text{N}^+\left(0, (A_{\tau}^{(m)})^2\right) \end{align}\]

\(A_0^{(m)}\) can be 0, implying that \(\beta_{u,1}^{(m)}\) is fixed at 0.

When \(U_m > 1\), sum-to-zero constraints (Section 5.1.2) can be applied.

5.4.2 Contribution to posterior density

\[\begin{equation} \text{N}(\tau_m \mid 0, A_{\tau}^{(m)2}) \prod_{u=1}^{U_m} \text{N}\left(\beta_{u,1}^{(m)} \mid 0, (A_0^{(m)})^2\right) \prod_{v=2}^{V_m} \text{N}\left(\beta_{u,v}^{(m)} \mid \beta_{u,v-1}^{(m)}, \tau_m^2 \right) \end{equation}\]

5.4.3 Forecasting

\[\begin{equation} \beta_{u,V_m+h}^{(m)} \sim \text{N}(\beta_{u,V_m+h-1}^{(m)}, \tau_m^2) \end{equation}\]

If the prior includes sum-to-zero constraints, means are subtracted from the forecasted values within each combination of ‘along’ and ‘by’ variables.

5.4.4 Code

RW(s = 1, along = NULL, zero_sum = TRUE)

• s is \(A_{\tau}^{(m)}\). Defaults to 1.
• sd is \(A_0^{(m)}\). Defaults to 1.
• along used to identify ‘along’ and ‘by’ dimensions.
• if zero_sum is TRUE, sum-to-zero constraints are applied.

5.5.1 Model

Second-order random walk

\[\begin{align} \beta_{u,1}^{(m)} & \sim \text{N}\left(0, (A_0^{(m)})^2\right) \\ \beta_{u,2}^{(m)} & \sim \text{N}\left(\beta_{u,1}, (A_{\eta}^{(m)})^2\right) \\ \beta_{u,v}^{(m)} & \sim \text{N}\left(2 \beta_{u,v-1}^{(m)} - \beta_{u,v-2}^{(m)}, \tau_m^2\right), \quad v = 3, \cdots, V_m \\ \tau_m & \sim \text{N}^+\left(0, (A_{\tau}^{(m)})^2\right) \end{align}\]

\(A_0^{(m)}\) can be 0, implying that \(\beta_{u,1}^{(m)}\) is fixed at 0.

When \(U_m > 1\), sum-to-zero constraints (Section 5.1.2) can be applied.

5.5.2 Contribution to posterior density

\[\begin{equation} \text{N}(\tau_m \mid 0, A_{\tau}^{(m)2}) \prod_{u=1}^{U_m} \text{N}(\beta_{u,1}^{(m)} \mid 0, (A_0^{(m)})^2) \text{N}(\beta_{u,2}^{(m)} \mid \beta_{u,1}^{(m)}, (A_{\eta}^{(m)})^2) \prod_{v=3}^{V_m} \text{N}\left(\beta_{u,v}^{(m)} \mid 2 \beta_{u,v-1}^{(m)} - \beta_{u,v-2}^{(m)}, \tau_m^2 \right) \end{equation}\]

5.5.3 Forecasting

\[\begin{equation} \beta_{u,V_m+h}^{(m)} \sim \text{N}(2 \beta_{u,V_m+h-1}^{(m)} - \beta_{u,V_m+h-2}^{(m)}, \tau_m^2) \end{equation}\]

If the prior includes sum-to-zero constraints, means are subtracted from the forecasted values within each combination of ‘along’ and ‘by’ variables.

5.5.4 Code

RW2(s = 1, sd = 1, sd_slope = 1, along = NULL, zero_sum = TRUE)

• s is \(A_{\tau}^{(m)}\)
• sd is \(A_0^{(m)}\)
• sd_slope is \(A_{\eta}^{(m)}\)
• along used to identify ‘along’ and ‘by’ dimensions
• if zero_sum is TRUE, sum-to-zero constraints are applied

5.6.1 Model

Random walk with seasonal effect

\[\begin{align} \beta_{u,v}^{(m)} & = \alpha_{u,v} + \lambda_{u,v}, \quad v = 1, \cdots, V_m \\ \alpha_{u,1}^{(m)} & \sim \text{N}\left(0, (A_0^{(m)})^2 \right) \\ \alpha_{u,v}^{(m)} & \sim \text{N}(\alpha_{u,v-1}^{(m)}, \tau_m^2), \quad v = 2, \cdots, V_m \\ \lambda_{u,v}^{(m)} & \sim \text{N}(0, A_{\lambda}^{(m)}), \quad v = 1, \cdots, S_m - 1 \\ \lambda_{u,v}^{(m)} & = -\sum_{s=1}^{S_m-1} \lambda_{u,v-s}^{(m)}, \quad v = S_m,\; 2S_m, \cdots \\ \lambda_{u,v}^{(m)} & \sim \text{N}(\lambda_{u,v-S_m}^{(m)}, \omega_m^2), \quad \text{otherwise} \\ \tau_m & \sim \text{N}^+\left(0, A_{\tau}^{(m)2}\right) \\ \omega_m & \sim \text{N}^+\left(0, A_{\omega}^{(m)2}\right) \end{align}\]

\(A_0^{(m)}\) can be 0, implying that \(\alpha_{u,1}^{(m)}\) is fixed at 0.

\(A_{\omega}^{(m)2}\) can be set to zero, implying that seasonal effects are fixed over time.

When \(U_m > 1\), sum-to-zero constraints (Section 5.1.2) can be applied.

5.6.2 Contribution to posterior density

\[\begin{equation} \begin{split} & \text{N}(\tau_m \mid 0, A_{\tau}^{(m)2}) \text{N}(\omega_m \mid 0, A_{\omega}^{(m)2}) \\ \quad \times & \prod_{u=1}^{U_m} \left( \text{N}\left(\alpha_{u,1}^{(m)} \mid 0, (A_0^{(m)})^2 \right) \prod_{v=2}^{V_m} \text{N}\left(\alpha_{u,v}^{(m)} \mid \alpha_{u,v-1}^{(m)}, \tau_m^2 \right) \prod_{v=1}^{S_m-1} \text{N}\left(\lambda_{u,v}^{(m)} \mid 0, (A_{\lambda}^{(m)})^2\right) \\ \prod_{v > S_m, \; (v-1) \bmod S_m \neq 0}^{V_m} \text{N}\left(\lambda_{u,v}^{(m)} \mid \lambda_{u,v-S_m}^{(m)}, \omega_m^2\right) \right) \end{split} \end{equation}\]

5.6.3 Forecasting

\[\begin{align} \alpha_{J_m+h}^{(m)} & \sim \text{N}(\alpha_{J_m+h-1}^{(m)}, \tau_m^2) \\ \lambda_{J_m+h}^{(m)} & \sim \text{N}(\lambda_{J_m+h-S_m}^{(m)}, \omega_m^2) \\ \beta_{J_m+h}^{(m)} & = \alpha_{J_m+h}^{(m)} + \lambda_{J_m+h}^{(m)} \end{align}\]

5.6.4 Code

RW_Seas(n_seas, s = 1, sd = 1, s_seas = 0, sd_seas = 1, along = NULL, zero_sum = FALSE)

• n_seas is \(S_m\)
• s is \(A_{\tau}^{(m)}\)
• sd is \(A_0^{(m)}\)
• s_seas is \(A_{\omega}^{(m)}\)
• sd_seas is \(A_{\lambda}^{(m)}\)
• along used to identify ‘along’ and ‘by’ dimensions
• if zero_sum is TRUE, sum-to-zero constraints are applied

5.7.1 Model

Second-order random work, with seasonal effect

\[\begin{align} \beta_{u,v}^{(m)} & = \alpha_{u,v} + \lambda_{u,v}, \quad v = 1, \cdots, V_m \\ \alpha_{u,1}^{(m)} & \sim \text{N}\left(0, (A_0^{(m)})^2 \right) \\ \alpha_{u,2}^{(m)} & \sim \text{N}\left(\alpha_{u,1}^{(m)}, (A_{\eta}^{(m)})^2\right) \\ \alpha_{u,v}^{(m)} & \sim \text{N}(2 \alpha_{u,v-1}^{(m)} - \alpha_{u,v-2}^{(m)}, \tau_m^2), \quad v = 3, \cdots, V_m \\ \lambda_{u,v}^{(m)} & \sim \text{N}(0, A_{\lambda}^{(m)}), \quad v = 1, \cdots, S_m - 1 \\ \lambda_{u,v}^{(m)} & = -\sum_{s=1}^{S_m-1} \lambda_{u,v}^{(m)}, \quad v = S_m,\; 2S_m, \cdots \\ \lambda_{u,v}^{(m)} & \sim \text{N}(\lambda_{u,v-S_m}^{(m)}, \omega_m^2), \quad \text{otherwise} \\ \tau_m & \sim \text{N}^+\left(0, A_{\tau}^{(m)2}\right) \\ \omega_m & \sim \text{N}^+\left(0, A_{\omega}^{(m)2}\right) \end{align}\]

\(A_0^{(m)}\) can be 0, implying that \(\alpha_{u,1}^{(m)}\) is fixed at 0.

\(A_{\omega}^{(m)2}\) can be set to zero, implying that seasonal effects are fixed over time.

When \(U_m > 1\), sum-to-zero constraints (Section 5.1.2) can be applied.

5.7.2 Contribution to posterior density

\[\begin{equation} \begin{split} & \text{N}(\tau_m \mid 0, A_{\tau}^{(m)2}) \text{N}(\omega_m \mid 0, A_{\omega}^{(m)2}) \\ \quad \times & \prod_{u=1}^{U_m} \left( \text{N}(\alpha_{u,1}^{(m)} \mid 0, (A_0^{(m)})^2 ) \text{N}(\alpha_{u,2}^{(m)} \mid \alpha_{u,1}^{(m)}, (A_{\eta}^{(m)})^2 ) \prod_{v=3}^{V_m} \text{N}(\alpha_{u,v}^{(m)} \mid 2 \alpha_{u,v-1}^{(m)} - \alpha_{u,v-2}^{(m)}, \tau_m^2 ) \\ \prod_{v=1}^{S_m-1} \text{N}(\lambda_{u,v}^{(m)} \mid 0, (A_{\lambda}^{(m)})^2) \prod_{v > S_m, (v-1) \bmod S_m \neq 0}^{V_m} \text{N}(\lambda_{u,v}^{(m)} \mid \lambda_{u,v-S_m}^{(m)}, \omega_m^2) \right) \end{split} \end{equation}\]

5.7.3 Forecasting

\[\begin{align} \alpha_{J_m+h}^{(m)} & \sim \text{N}(2 \alpha_{J_m+h-1}^{(m)} - \alpha_{J_m+h-2}^{(m)}, \tau_m^2) \\ \lambda_{J_m+h}^{(m)} & \sim \text{N}(\lambda_{J_m+h-S_m}^{(m)}, \omega_m^2) \\ \beta_{J_m+h}^{(m)} & = \alpha_{J_m+h}^{(m)} + \lambda_{J_m+h}^{(m)} \end{align}\]

5.7.4 Code

RW2_Seas(n_seas, s = 1, sd = 1, sd_slope = 1, s_seas = 0, sd_seas = 1, along = NULL, zero_sum = FALSE)

• n_seas is \(S_m\)
• s is \(A_{\tau}^{(m)}\)
• sd is \(A_0^{(m)}\)
• sd_slope is \(A_{\eta}^{(m)}\)
• s_seas is \(A_{\omega}^{(m)}\)
• sd_seas is \(A_{\lambda}^{(m)}\)
• along used to identify ‘along’ and ‘by’ dimensions
• if zero_sum is TRUE, sum-to-zero constraints are applied

5.8.1 Model

\[\begin{equation} \beta_{u,v}^{(m)} \sim \text{N}\left(\phi_1^{(m)} \beta_{u,v-1}^{(m)} + \cdots + \phi_{K_m}^{(m)} \beta_{u,v-{K_m}}^{(m)}, \omega_m^2\right), \quad v = K_m + 1, \cdots, V_m. \end{equation}\] Internally, TMB derives values for \(\beta_{u,v}^{(m)}, v = 1, \cdots, K_m\), and for \(\omega_m\), that imply a stationary distribution, and that give every term \(\beta_{u,v}^{(m)}\) the same marginal variance. We denote this marginal variance \(\tau_m^2\), and assign it a prior \[\begin{equation} \tau_m \sim \text{N}^+(0, A_{\tau}^{(m)2}). \end{equation}\] Each of the \(\phi_k^{(m)}\) has prior \[\begin{equation} \frac{\phi_k^{(m)} + 1}{2} \sim \text{Beta}(S_1^{(m)}, S_2^{(m)}). \end{equation}\]

5.8.2 Contribution to posterior density

\[\begin{equation} \text{N}^+\left(\tau_m \mid 0, A_{\tau}^{(m)2} \right) \prod_{k=1}^{K_m} \text{Beta}\left(\tfrac{1}{2} \phi_k^{(m)} + \tfrac{1}{2} \mid 2, 2 \right) \prod_{u=1}^{U_m} p\left( \beta_{u,1}^{(m)}, \cdots, \beta_{u,V_m}^{(m)} \mid \phi_1^{(m)}, \cdots, \phi_{K_m}^{(m)}, \tau_m \right) \end{equation}\] where \(p\left( \beta_{u,1}^{(m)}, \cdots, \beta_{u,V_m}^{(m)} \mid \phi_1^{(m)}, \cdots, \phi_{K_m}^{(m)}, \tau_m \right)\) is calculated internally by TMB.

5.8.3 Forecasting

\[\begin{equation} \beta_{u,V_m + h}^{(m)} \sim \text{N}\left(\phi_1^{(m)} \beta_{u,V_m + h - 1}^{(m)} + \cdots + \phi_{K_m}^{(m)} \beta_{u,V_m+h-K_m}^{(m)}, \tau_m^2\right) \end{equation}\]

5.8.4 Code

AR(n_coef = 2,
   s = 1,
   shape1 = 5,
   shape2 = 5,
   along = NULL,
   zero_sum = FALSE)

• n_coef is \(K_m\)
• s is \(A_{\tau}^{(m)}\)
• shape1 is \(S_1^{(m)}\)
• shape2 is \(S_2^{(m)}\)
• along is used to indentify the ‘along’ and ‘by’ dimensions

5.9.1 Model

\[\begin{align} \beta_{u,1}^{(m)} & \sim \text{N}(0, \tau_m^2), \quad u = 1, \cdots, U_m \\ \beta_{u,v}^{(m)} & \sim \text{N}(\phi_m \beta_{u,v-1}^{(m)}, (1 - \phi_m^2) \tau_m^2), \quad u = 1, \cdots, U_m, \quad v = 2, \cdots, V_m \\ \phi_m & = a_{0,m} + (a_{1,m} - a_{0,m}) \phi_m^{\prime} \\ \phi_m^{\prime} & \sim \text{Beta}(S_1^{(m)}, S_2^{(m)}) \\ \tau_m & \sim \text{N}^+\left(0, A_{\tau}^{(m)2}\right). \end{align}\] This is adapted from the specification used for AR1 densities in TMB. It implies that the marginal variance of all \(\beta_{u,v}^{(m)}\) is \(\tau_m^2\). We require that \(-1 < a_{0m} < a_{1m} < 1\).

5.9.2 Contribution to posterior density

\[\begin{equation} \text{N}(\tau_m \mid 0, A_{\tau}^{(m)2}) \text{Beta}( \phi_m^{\prime} \mid S_1^{(m)}, S_2^{(m)}) \prod_{u=1}^{U_m} \text{N}\left(\beta_{u,1}^{(m)} \mid 0, \tau_m^2 \right) \prod_{u=1}^{U_m} \prod_{j=2}^{V_m} \text{N}\left(\beta_{u,v}^{(m)} \mid \phi_m \beta_{u,v-1}^{(m)}, (1 - \phi_m^2) \tau_m^2 \right) \end{equation}\]

5.9.3 Forecasting

\[\begin{equation} \beta_{J_m + h}^{(m)} \sim \text{N}\left(\phi_m \beta_{J_m + h - 1}^{(m)}, (1 - \phi_m^2) \tau_m^2\right) \end{equation}\]

5.9.4 Code

AR1(s = 1,
    shape1 = 5,
    shape2 = 5,
    min = 0.8,
    max = 0.98,
    along = NULL,
    zero_sum = TRUE)

• s is \(A_{\tau}^{(m)}\)
• shape1 is \(S_1^{(m)}\)
• shape2 is \(S_2^{(m)}\)
• min is \(a_{0m}\)
• max is \(a_{1m}\)
• along is used to identify ‘along’ and ‘by’ dimensions

The defaults for min and max are based on the defaults for function ets() in R package forecast (Hyndman and Khandakar 2008).

5.10.1 Model

\[\begin{align} \beta_{u,v}^{(m)} & = \alpha_{u,v}^{(m)} + \epsilon_{u,v}^{(m)} \\ \alpha_{u,v}^{(m)} & = \left(v - \frac{V_m + 1}{2}\right) \eta_u^{(m)} \\ \eta_u^{(m)} & \sim \text{N}\left(B_{\eta}^{(m)}, (A_{\eta}^{(m)})^2\right)\\ \epsilon_{u,v}^{(m)} & \sim \text{N}(0, \tau_m^2) \\ \tau_m & \sim \text{N}^+\left(0, (A_{\tau}^{(m)})^2\right) \end{align}\]

Note that \(\sum_{v=1}^{V_m} \alpha_{u,v}^{(m)} = 0\).

5.10.2 Contribution to posterior density

\[\begin{equation} \text{N}(\tau_m \mid 0, A_{\tau}^{(m)2}) \prod_{u=1}^{U_m} \text{N}(\eta_u^{(m)} \mid B_{\eta}^{(m)}, A_{\eta}^{(m)2}) \prod_{u=1}^{U_m} \prod_{v=1}^{V_m} \text{N}\left(\beta_{u,v}^{(m)} \mid v - \frac{V_m + 1}{2}, \tau_m^2 \right) \end{equation}\]

5.10.3 Forecasting

\[\begin{equation} \beta_{u,V_m + h}^{(m)} \sim \text{N}\left(\left(\frac{V_m - 1}{2}+ h\right) \eta_u^{(m)}, \tau_m^2\right) \end{equation}\]

5.10.4 Code

Lin(s = 1, mean_slop = 0, sd_slope = 1, along = NULL, zero_sum = FALSE)

• s is \(A_{\tau}^{(m)}\)
• mean_slope is \(B_{\eta}^{(m)}\)
• sd_slope is \(A_{\eta}^{(m)}\)
• along is used to indentify ‘along’ and ‘by’ dimensions

5.11.1 Model

\[\begin{align} \beta_{u,v}^{(m)} & = \alpha_{u,v}^{(m)} + \epsilon_{u,v}^{(m)} \\ \alpha_{u,v}^{(m)} & = \left(v - \frac{V_m + 1}{2}\right) \eta_u^{(m)} \\ \eta_u^{(m)} & \sim \text{N}\left(B_{\eta}^{(m)}, (A_{\eta}^{(m)})^2\right)\\ \epsilon_{u,v}^{(m)} & \sim \text{N}\left(\phi_1^{(m)} \epsilon_{u,v-1}^{(m)} + \cdots + \phi_{K_m}^{(m)} \epsilon_{u,v-{K_m}}^{(m)}, \omega_m^2\right), \quad v = K_m + 1, \cdots, V_m \end{align}\]

Note that \(\sum_{v=1}^{V_m} \alpha_{u,v}^{(m)} = 0\).

Internally, TMB derives values for \(\epsilon_{u,v}^{(m)}, v = 1, \cdots, K_m\), and for \(\omega_m\), that provide the \(\epsilon_{u,v}^{(m)}\) with a stationary distribution in which each term has the same marginal variance. We denote this marginal variance \(\tau_m^2\), and assign it a prior \[\begin{equation} \tau_m \sim \text{N}^+(0, A_{\tau}^{(m)2}). \end{equation}\] Each of the individual \(\phi_k^{(m)}\) has prior \[\begin{equation} \frac{\phi_k^{(m)} + 1}{2} \sim \text{Beta}(S_1^{(m)}, S_2^{(m)}). \end{equation}\]

5.11.2 Contribution to posterior density

\[\begin{equation} \begin{split} & \text{N}^+\left(\tau_m \mid 0, A_{\tau}^{(m)2} \right) \prod_{k=1}^{K_m} \text{Beta}\left( \tfrac{1}{2} \phi_k^{(m)} + \tfrac{1}{2} \mid 2, 2 \right) \\ & \quad \times \prod_{u=1}^{U_m} \text{N}(\eta_u^{(m)} \mid 0, A_{\eta}^{(m)2}) p\left( \epsilon_{u,1}^{(m)}, \cdots, \epsilon_{u,V_m}^{(m)} \mid \phi_1^{(m)}, \cdots, \phi_{K_m}^{(m)}, \tau_m \right) \end{split} \end{equation}\] where \(p\left( \epsilon_{u,1}^{(m)}, \cdots, \epsilon_{u,V_m}^{(m)} \mid \phi_1^{(m)}, \cdots, \phi_{K_m}^{(m)}, \tau_m \right)\) is calculated internally by TMB.

5.11.3 Forecasting

\[\begin{align} \beta_{u, V_m + h}^{(m)} & = \left(\frac{V_m - 1}{2}+ h\right) \eta_u^{(m)} + \epsilon_{u,V_m+h}^{(m)} \\ \epsilon_{u,V_m+h}^{(m)} & \sim \text{N}\left(\phi_1^{(m)} \epsilon_{u,V_m + h - 1}^{(m)} + \cdots + \phi_{K_m}^{(m)} \epsilon_{u,V_m+h-K_m}^{(m)}, \omega_m^2\right) \end{align}\]

5.11.4 Code

Lin_AR(n_coef = 2,
       s = 1,
       shape1 = 5,
       shape2 = 5,
       mean_slope = 0,
       sd_slope = 1,
       along = NULL,
       zero_coef = FALSE)

• n_coef is \(K_m\)
• s is \(A_{\tau}^{(m)}\)
• shape1 is \(S_1^{(m)}\)
• shape2 is \(S_2^{(m)}\)
• mean_slope is \(B_{\eta}^{(m)}\)
• sd_slope is \(A_{\eta}^{(m)}\)
• along is used to indentify ‘along’ and ‘by’ variables

5.13.1 Model

Penalised spline (P-spline)

\[\begin{equation} \pmb{\beta}_u^{(m)} = \pmb{B}^{(m)} \pmb{\alpha}_u^{(m)}, \quad u = 1, \cdots, U_m \end{equation}\] where \(\pmb{\beta}_u^{(m)}\) is the subvector of \(\pmb{\beta}^{(m)}\) composed of elements from the \(u\)th combination of the ‘by’ variables, \(\pmb{B}^{(m)}\) is a \(V_m \times K_m\) matrix of B-splines, and \(\pmb{\alpha}_u^{(m)}\) has a second-order random walk prior (Section 5.5).

\(\pmb{B}^{(m)} = (\pmb{b}_1^{(m)}(\pmb{v}), \cdots, \pmb{b}_{K_m}^{(m)}(\pmb{v}))\), with \(\pmb{v} = (1, \cdots, V_m)^{\top}\). The B-splines are centered, so that \(\pmb{1}^{\top} \pmb{b}_k^{(m)}(\pmb{v}) = 0\), \(k = 1, \cdots, K_m\).

5.13.2 Contribution to posterior density

\[\begin{equation} \text{N}(\tau_m \mid 0, A_{\tau}^{(m)2}) \prod_{u=1}^{U_m} \prod_{k=1}^2 \text{N}(\alpha_{u,k}^{(m)} \mid 0, 1) \prod_{u=1}^{U_m}\prod_{k=3}^{K_m} \text{N}\left(\alpha_{u,k}^{(m)} - 2 \alpha_{u,k-1}^{(m)} + \alpha_{u,k-2}^{(m)} \mid 0, \tau_m^2 \right) \end{equation}\]

5.13.3 Forecasting

Terms with a P-Spline prior cannot be forecasted.

5.13.4 Code

Sp(n = NULL, s = 1)

• n is \(K_m\). Defaults to \(\max(0.7 J_m, 4)\).
• s is the \(A_{\tau}^{(m)}\) from the second-order random walk prior. Defaults to 1.
• along is used to identify ‘along’ and ‘by’ variables

5.14.1 Model

Age but no sex or gender

Let \(\pmb{\beta}_u\) be the age effect for the \(u\)th combination of the ‘by’ variables. With an SVD prior, \[\begin{equation} \pmb{\beta}_u^{(m)} = \pmb{F}^{(m)} \pmb{\alpha}_u^{(m)} + \pmb{g}^{(m)}, \quad u = 1, \cdots, U_m \end{equation}\] where \(\pmb{F}^{(m)}\) is a \(V_m \times K_m\) matrix, and \(\pmb{g}^{(m)}\) is a vector with \(V_m\) elements, both derived from a singular value decomposition (SVD) of an external dataset of age-specific values for all sexes/genders combined. The construction of \(\pmb{F}^{(m)}\) and \(\pmb{g}^{(m)}\) is described in Appendix 13.2. The centering and scaling used in the construction allow use of the simple prior \[\begin{equation} \alpha_{u,k}^{(m)} \sim \text{N}(0, 1), \quad u = 1, \cdots, U_m, k = 1, \cdots, K_m. \end{equation}\]

Joint model of age and sex/gender

In the joint model, vector \(\pmb{\beta}_u\) represents the interaction between age and sex/gender for the \(u\)th combination of the ‘by’ variables. Matrix \(\pmb{F}^{(m)}\) and vector \(\pmb{g}^{(m)}\) are calculated from data that separate sexes/genders. The model is otherwise unchanged.

Independent models for each sex/gender

In the independent model, vector \(\pmb{\beta}_{s,u}\) represents age effects for sex/gender \(s\) and the \(u\)th combination of the ‘by’ variables, and we have \[\begin{equation} \pmb{\beta}_{s,u}^{(m)} = \pmb{F}_s^{(m)} \pmb{\alpha}_{s,u}^{(m)} + \pmb{g}_s^{(m)}, \quad s = 1, \cdots, S; \quad u = 1, \cdots, U_m \end{equation}\] Matrix \(\pmb{F}_s^{(m)}\) and vector \(\pmb{g}_s^{(m)}\) are calculated from data that separate sexes/genders. The prior is \[\begin{equation} \alpha_{s,u,k}^{(m)} \sim \text{N}(0, 1), \quad s = 1, \cdots, S; \quad u = 1, \cdots, U_m; \quad k = 1, \cdots, K_m. \end{equation}\]

5.14.2 Contribution to posterior density

\[\begin{equation} \prod_{u=1}^{U_m}\prod_{k=1}^{K_m} \text{N}\left(\alpha_{uk}^{(m)} \mid 0, 1 \right) \end{equation}\] for the age-only and joint models, and \[\begin{equation} \prod_{s=1}^S \prod_{u=1}^{U_m}\prod_{k=1}^{K_m} \text{N}\left(\alpha_{s,u,k}^{(m)} \mid 0, 1 \right) \end{equation}\] for the independent model

5.14.3 Forecasting

Terms with an SVD prior cannot be forecasted.

5.14.4 Code

SVD(ssvd, n_comp = NULL, indep = TRUE)

where - ssvd is an object containing \(\pmb{F}\) and \(\pmb{g}\) - n_comp is the number of components to be used (which defaults to ceiling(n/2), where n is the number of components in ssvd - indep determines whether and independent or joint model will be used if the term being modelled contains a sex or gender variable.

5.15.1 Model

The SVD_RW() prior is identical to the SVD() prior except that the coefficients evolve over time, following independent random walks. For instance, in the combined-sex/gender and joint models with \(K_m\) SVD components,

\[\begin{align} \pmb{\beta}_{u,t}^{(m)} & = \pmb{F}^{(m)} \pmb{\alpha}_{u,t}^{(m)} + \pmb{g}^{(m)} \\ \alpha_{u,k,1}^{(m)} & \sim \text{N}(0, (A_0^{(m)})^2), \\ \alpha_{u,k,t}^{(m)} & \sim \text{N}(\alpha_{u,k,t-1}^{(m)}, \tau_m^2), \quad t = 2, \cdots, T \\ \tau_m & \sim \text{N}^+\left(0, (A_{\tau}^{(m)})^2\right) \end{align}\]

5.15.2 Contribution to posterior density

In the combined-sex/gender and joint models,

\[\begin{equation} \text{N}(\tau_m \mid 0, A_{\tau}^{(m)2}) \prod_{u=1}^{U_m} \prod_{k=1}^{K_m} \text{N}(\alpha_{u,k,1}^{(m)} \mid 0, (A_0^{(m)})^2) \prod_{t=2}^{T} \text{N}\left(\alpha_{u,k,t}^{(m)} \mid \alpha_{u,k,t-1}^{(m)}, \tau_m^2 \right), \end{equation}\]

and in the independent model,

\[\begin{equation} \text{N}(\tau_m \mid 0, A_{\tau}^{(m)2}) \prod_{u=1}^{U_m} \prod_{s=1}^{S} \prod_{k=1}^{K_m} \text{N}(\alpha_{u,s,k,1}^{(m)} \mid 0, (A_0^{(m)})^2) \prod_{t=2}^{T} \text{N}\left(\alpha_{u,s,k,t}^{(m)} \mid \alpha_{u,s,k,t-1}^{(m)}, \tau_m^2 \right) \end{equation}\]

5.15.3 Forecasting

\[\begin{align} \alpha_{u,k,T+h}^{(m)} & \sim \text{N}(\alpha_{u,k,T+h-1}^{(m)}, \tau_m^2) \\ \pmb{\beta}_{u,T+h}^{(m)} & = \pmb{F}^{(m)} \pmb{\alpha}_{u,T+h}^{(m)} + \pmb{g}^{(m)} \end{align}\]

5.15.4 Code

SVD_RW(ssvd,
       n_comp = NULL,
       indep = TRUE,
       s = 1,
       sd = 1,
       zero_sum = FALSE)

where - ssvd is an object containing \(\pmb{F}\) and \(\pmb{g}\) - n_comp is \(K_m\) - indep determines whether and independent or joint model will be used if the term being modelled contains a sex or gender variable. - s is \(A_{\tau}^{(m)}\) - sd is \(A_0^{(m)}\)

5.19.1 Model

Elements of \(\pmb{\beta}^{(m)}\) are treated as known with certainty.

5.19.2 Contribution to posterior density

Known priors make no contribution to the posterior density.

5.19.3 Forecasting

Main effects with a known prior cannot be forecasted.

5.19.4 Code

Known(values)

• values is a vector containing the \(\beta_j^{(m)}\).

6.1.1 Standard prior

\[\begin{align} \zeta_p \mid \varphi & \sim \text{N}(0, \varphi^2) \\ \varphi & \sim \text{N}^+(0, 1) \end{align}\]

6.1.2 Shrinkage prior

Regularized horseshoe prior (Piironen and Vehtari 2017)

\[\begin{align} \zeta_p \mid \vartheta_p, \varphi & \sim \text{N}(0, \vartheta_p^2 \varphi^2) \\ \vartheta_p & \sim \text{Cauchy}^+(0, 1) \\ \varphi & \sim \text{Cauchy}^+(0, A_{\varphi}^2) \\ A_{\varphi} & = \frac{p_0}{p_0 + P} \frac{\hat{\sigma}}{\sqrt{n}} \end{align}\] where \(p_0\) is an initial guess at the number of \(\zeta_p\) that are non-zero, and \(\hat{\sigma}\) is obtained as follows:

• Poisson. Using maximum likelihood, fit the GLM \[\begin{align} y_i & \sim \text{Poisson}(w_i \gamma_i) \\ \log \gamma_i & = \sum_{m=0}^{M} \pmb{X}^{(m)} \pmb{\beta}^{(m)}, \end{align}\] and set \[\begin{equation} \hat{\sigma} = \frac{1}{n}\sum_{i=1}^n \frac{1}{w_i \hat{\gamma}_i}. \end{equation}\]
• Binomial. Using maximum likelihood, fit the GLM \[\begin{align} y_i & \sim \text{binomial}(w_i, \gamma_i) \\ \text{logit} \gamma_i & = \sum_{m=0}^{M} \pmb{X}^{(m)} \pmb{\beta}^{(m)}, \end{align}\] and set \[\begin{equation} \hat{\sigma} = \frac{1}{n}\sum_{i=1}^n \frac{1}{w_i \hat{\gamma}_i (1 - \hat{\gamma}_i)}. \end{equation}\]
• Normal. Using maximum likelihood, fit the linear model \[\begin{equation} y_i \sim \text{N}\left(\sum_{m=0}^{M} \pmb{X}^{(m)} \pmb{\beta}^{(m)}, w_i^{-1}\xi^2 \right) \end{equation}\] and set \(\hat{\sigma} = \hat{\xi}\).

The quantities used for Poisson and binomial likelihoods are derived from normal approximations to GLMs (Piironen and Vehtari 2017; Gelman et al. 2014, sec. 16.2).

8.1.1 Random Rounding to Base 3

Random rounding to base 3 (RR3) is a confidentialization method used by some statistical agencies. It is applied to counts data. Each count \(x\) is rounded randomly as follows:

• If \(x \mod 3 = 0\), then \(x\) is left unchanged;
• if \(x \mod 3 = 1\) then \(x\) is changed to \(x-1\) with probability 2/3, and is changed to \(x + 2\) with probability 1/3; and
• if \(x \mod 3 = 2\) then \(x\) is changed to \(x-2\) with probability 1/3, and is changed to \(x + 1\) with probability 2/3.

RR3 data models can be used with Poisson or binomial likelihoods. Let \(y_i\) denote the observed value for the outcome, and \(y_i^*\) the true value. The likelihood with a RR3 data model is then

\[\begin{align} p(y_i | \gamma_i, w_i) & = \sum_{y_i^*} p(y_i | y_i^*) p(y_i^* | \gamma_i, w_i) \\ & = \sum_{k = -2, -1, 0, 1, 2} p(y_i | y_i + k) p(y_i + k | \gamma_i, w_i) \\ & = \tfrac{1}{3} p(y_i - 2 | \gamma_i, w_i) + \tfrac{2}{3} p(y_i - 1 | \gamma_i, w_i) + p(y_i | \gamma_i, w_i) + \tfrac{2}{3} p(y_i + 1 | \gamma_i, w_i) + \tfrac{1}{3} p(y_i + 2 | \gamma_i, w_i) \end{align}\]

9.2.1 Step 0: Select ‘inner’ and ‘outer’ variables

Select variables to be used in inner model. By default, these are the age, sex, and time variables in the model. All remaining variables are ‘outer’ variables.

9.2.2 Step 1: Fit inner model

Aggregate the data using the classification formed by the inner variables. (See Section 9.1 on aggregation procedures.) Remove all terms not involving ‘inner’ variables, other than the intercept term, from the model. Set dispersion to 0. Fit the resulting model.

9.2.3 Step 2: Fit outer model

Let \(\hat{\mu}_i^{\text{in}}\) be point estimates for the linear predictor \(\mu_i\) obtained from the inner model.

Poisson model

Aggregate the data using the classification formed by the outer variables. Remove all terms involving the ‘inner’ variables, plus the intercept, from the model. Set dispersion to 0. Set exposure to \(w_i^{\text{out}} = \hat{\mu}^{\text{in}} w_i\). Fit the model.

Binomial model

Fit the original model, but set dispersion to 0, and for all terms from the ‘inner’ model, use Known priors using point estimates from the inner model.

Normal model

Aggregate the data using the classification formed by the outer variables. Remove all terms involving the ‘inner’ variables, plus the intercept, from the model. Set the outcome variable to to $y_i^{} = y_i - ^{}. Fit the model.

9.2.4 Step 4: Concatenate estimates

Concatenate posterior distributions for the inner terms from the inner model to posterior distributions for the outer terms from the outer model.

9.2.5 Step 5: Calculate dispersion

If the original model includes a dispersion term, then estimate dispersion. Let \(\hat{\mu}_i^{\text{comb}}\) be point estimates for the linear predictor obtained from the concatenated estimates.

Poisson model

Use the original disaggregated data, or, if the original data contains more then 10,000 rows, select 10,000 rows at random from the original data. Remove all terms from the original model except for the intercept. Set exposure to \(w_i^{\text{out}} = \hat{\mu}^{\text{comb}} w_i\).

Binomial model

Fit the the original model, but with all terms except the intercept having Known priors, where the values are obtained from point estimates from the concatenated estimates.

Normal model

Use the original disaggregated data, or, if the original data contains more then 10,000 rows, select 10,000 rows at random from the original data. Remove all terms from the original model except for the intercept. Set the outcome to \(y_i^{\text{out}} = y_i - \hat{\mu}^{\text{comb}}\).

12.1.1 Poisson likelihood

12.1.2 Binomial likelihood

12.1.3 Normal likelihood

\[\begin{equation} y_i^{\text{rep}} \sim \text{N}(\gamma_i, \xi^2 / w_i) \end{equation}\]

12.1.4 Data models for outcomes

If the overall model includes a data model for the outcome, then a further set of draws is made, deriving values for the observed outcomes, given values for the true outcomes.