library(bage)Specification document - a mathematical description of models used by bage.
Note: some features described here have not been implemented yet.
Let \(y_i\) be a count of events in cell \(i = 1, \cdots, n\) and let \(w_i\) be the corresponding exposure measure, with the possibility that \(w_i \equiv 1\). The likelihood under the Poisson model is then \[\begin{align} y_i & \sim \text{Poisson}(\gamma_i w_i) \tag{2.1} \\ \gamma_i & \sim \text{Gamma}\left(\xi^{-1}, (\mu_i \xi)^{-1}\right), \tag{2.2} \end{align}\] using the shape-rates parameterisation of the Gamma distribution. Parameter \(\xi\) governs dispersion, with \[\begin{equation} \text{var}(\gamma_i \mid \mu_i, \xi) = \xi \mu_i^2 \end{equation}\] and \[\begin{equation} \text{var}(y_i \mid \mu_i, \xi, w_i) = (1 + \xi \mu_i w_i ) \times \mu_i w_i. \end{equation}\] We allow \(\xi\) to equal 0, in which case the model reduces to \[\begin{equation} y_i \sim \text{Poisson}(\mu_i w_i). \end{equation}\]
For \(\xi > 0\), Equations (2.1) and (2.2) are equivalent to \[\begin{equation} y_i \sim \text{NegBinom}\left(\xi^{-1}, (1 + \mu_i w_i \xi)^{-1}\right) \end{equation}\] (Norton, Christen, and Fox 2018; Simpson 2022). This is the format we use internally for estimation. When values for \(\gamma_i\) are needed, we generate them on the fly, using the fact that \[\begin{equation} \gamma_i \mid y_i, w_i, \mu_i, \xi \sim \text{Gamma}\left(y_i + \xi^{-1}, w_i + (\xi \mu_i)^{-1}\right). \end{equation}\]
The likelihood under the binomial model is \[\begin{align} y_i & \sim \text{Binomial}(w_i, \gamma_i) \tag{2.3} \\ \gamma_i & \sim \text{Beta}\left(\xi^{-1} \mu_i, \xi^{-1}(1 - \mu_i)\right). \tag{2.4} \end{align}\] Parameter \(\xi\) again governs dispersion, with \[\begin{equation} \text{var}(\gamma_i \mid \mu_i, \xi) = \frac{\xi}{1 + \xi} \times \mu_i (1 -\mu_i) \end{equation}\] and \[\begin{equation} \text{var}(y_i \mid w_i, \mu_i, \xi) = \frac{\xi w_i + 1}{\xi + 1} \times w_i \mu_i (1 - \mu_i). \end{equation}\]
We allow \(\xi\) to equal 0, in which case the model reduces to \[\begin{equation} y_i \text{Binom}(w_i, \mu_i). \end{equation}\] Equations (2.3) and (2.4) are equivalent to \[\begin{equation} y_i \sim \text{BetaBinom}\left(w_i, \xi^{-1} \mu_i, \xi^{-1} (1 - \mu_i) \right), \end{equation}\] which is what we use internally. Values for \(\gamma_i\) can be generated using \[\begin{equation} \gamma_i \mid y_i, w_i, \mu_i, \xi \sim \text{Beta}\left(y_i + \xi^{-1} \mu_i, w_i - y_i + \xi^{-1}(1-\mu_i) \right). \end{equation}\]
\[\begin{equation} y_i \sim \text{N}(\mu_i, w_i^{-1}\xi^2) \end{equation}\] where the \(w_i\) are weights.
Response \(y_i\) is standardized to have mean 0 and standard deviation 1. We set \[\begin{equation} y_i = \frac{y_i^{*} - \bar{y}^*}{s^*} \end{equation}\] where the \(y_i^*\) are the values originally supplied by the user, and \(\bar{y}^*\) and \(s^*\) are the mean and standard deviation of the \(y_i^*\).
Let \(\pmb{\mu} = (\mu_1, \cdots, \mu_n)^{\top}\). Our model for \(\pmb{\mu}\) is \[\begin{equation} \pmb{\mu} = \sum_{m=0}^{M} \pmb{X}^{(m)} \pmb{\beta}^{(m)} + \pmb{Z} \pmb{\zeta} \tag{3.1} \end{equation}\] where
The algorithm for assigning default priors:
The intercept term \(\pmb{\beta}^{(0)}\) can only be given a fixed-normal prior (Section 4.2) or a Known prior (Section 4.17).
\[\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}\]
\[\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}\]
\[\begin{align} \tau_m & \sim \text{N}^+\left(0, A_{\tau}^{(m)2}\right) \\ \beta_j^{(m)} & \sim \text{N}\left(0, \tau_m^2 \right) \end{align}\]
\[\begin{equation} \beta_{J_m+h+1}^{(m)} \sim \text{N}(0, \tau_m^2) \end{equation}\]
N(s = 1)s is \(A_{\tau}^{(m)}\). Defaults to 1.\[\begin{equation} \beta_j^{(m)} \sim \text{N}\left(0, A_{\beta}^{(m)2}\right) \end{equation}\]
\[\begin{equation} \prod_{j=1}^{J_m} \text{N}(\beta_j^{(m)} \mid 0, A_{\beta}^{(m)2}) \end{equation}\]
\[\begin{equation} \beta_j^{(m)} \sim \text{N}\left(0, A_{\beta}^{(m)2}\right) \end{equation}\]
\[\begin{equation} \beta_{J_m+h+1}^{(m)} \sim \text{N}(0, A_{\beta}^{(m)2}) \end{equation}\]
NFix(sd = 1)sd is \(A_{\tau}^{(m)}\). Defaults to 1.\[\begin{align} \beta_{u1}^{(m)} & \sim \text{N}(0, 1), \quad u = 1, \cdots, U_m \\ \beta_{uv}^{(m)} & \sim \text{N}(\beta_{u,v-1}^{(m)}, \tau_m^2), \quad u = 1, \cdots, U_m, v = 2, \cdots, V_m \\ \tau_m & \sim \text{N}^+\left(0, A_{\tau}^{(m)2}\right) \end{align}\]
\[\begin{equation} \text{N}(\tau_m \mid 0, A_{\tau}^{(m)2}) \prod_{u=1}^{U_m} \text{N}(\beta_{u1}^{(m)} \mid 0, 1) \prod_{u=1}^{U_m} \prod_{v=2}^{V_m} \text{N}\left(\beta_{uv}^{(m)} \mid \beta_{u,v-1}^{(m)}, \tau_m^2 \right) \end{equation}\]
\[\begin{align} \tau_m & \sim \text{N}^+\left(0, A_{\tau}^{(m)2}\right) \\ \beta_{u,1}^{(m)} & \sim \text{N}(0, 1), \quad u = 1, \cdots, U_m, \\ \beta_{u,v}^{(m)} & \sim \text{N}(\beta_{u,v-1}^{(m)}, \tau_m^2), \quad u = 1, \cdots, U_m, v = 2, \cdots, V_m \end{align}\]
\[\begin{equation} \beta_{u,V_m+h}^{(m)} \sim \text{N}(\beta_{u,V_m+h-1}^{(m)}, \tau_m^2), \quad u = 1, \cdots, U_m \end{equation}\]
RW(s = 1, along = NULL)s is \(A_{\tau}^{(m)}\). Defaults to 1.along used to identify “along” and “by” dimensions\[\begin{align} \beta_{u,v}^{(m)} & \sim \text{N}(0, 1), \quad u = 1,\cdots, U_m, \quad v = 1,2 \\ \beta_{u,v}^{(m)} & \sim \text{N}(2 \beta_{u,v-1}^{(m)} - \beta_{u,v-2}^{(m)}, \tau_m^2), \quad u = 1, \cdots, U_m, \quad v = 3, \cdots, V_m \\ \tau_m & \sim \text{N}^+\left(0, A_{\tau}^{(m)2}\right) \end{align}\]
\[\begin{equation} \text{N}(\tau_m \mid 0, A_{\tau}^{(m)2}) \prod_{u=1}^{U_m} \prod_{v=1}^2 \text{N}(\beta_{u,v}^{(m)} \mid 0, 1) \prod_{u=1}^{U_m}\prod_{v=3}^{V_m} \text{N}\left(\beta_{u,v}^{(m)} - 2 \beta_{u,v-1}^{(m)} + \beta_{u,v-2}^{(m)} \mid 0, \tau_m^2 \right) \end{equation}\]
\[\begin{align} \tau_m & \sim \text{N}^+\left(0, A_{\tau}^{(m)2}\right) \\ \beta_{u,v}^{(m)} & \sim \text{N}(0, 1), \quad u = 1, \cdots, U_m, \quad v = 1,2 \\ \beta_{u,v}^{(m)} & \sim \text{N}(2 \beta_{u,v-1}^{(m)} - \beta_{u,v-2}^{(m)}, \tau_m^2), \quad u = 1, \cdots, U_m, \quad v = 3,\cdots,V_m \end{align}\]
\[\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}\]
RW2(s = 1, sd = 1, along = NULL)s is \(A_{\tau}^{(m)}\)sd is \(A_{\eta}^{(m)}\)along used to identify “along” and “by” dimensions\[\begin{align} \beta_{u,v}^{(m)} & = \alpha_{u,v} + \lambda_{u,s_v}, \quad u = 1, \cdots, U_m, \quad v = 1, \cdots, V_m \\ \alpha_{u,1}^{(m)} & \sim \text{N}(0, 1), \quad u = 1, \cdots, U_m \\ \alpha_{u,v}^{(m)} & \sim \text{N}(\alpha_{u,v-1}^{(m)}, \tau_m^2), \quad u = 1, \cdots, U_m, \quad v = 2, \cdots, V_m \\ \lambda_{u,v}^{(m)} & \sim \text{N}(0, 1), \quad u = 1, \cdots, U_m, \quad v = 1, \cdots, S_m \\ \lambda_{u,v}^{(m)} & \sim \text{N}(\lambda_{u,v-S_m}^{(m)}, \omega_m^2), \quad u = 1, \cdots, U_m, \quad v = S_m + 1, \cdots, V_m \tag{4.1} \\ \\ \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}\]
We allow \(A_{\omega}^{(m)2}\) to be set to zero, in which case (4.1) reduces to \[\begin{equation} \lambda_{u,v}^{(m)} = \lambda_{u,v-S_m}^{(m)}, \quad u = 1, \cdots, U_m, \quad v = S_m + 1, \cdots, V_m, \end{equation}\] implying that seasonal effects are constant across years.
\[\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, 1 ) \prod_{v=2}^{V_m} \text{N}(\alpha_{u,v}^{(m)} \mid \alpha_{u,v-1}^{(m)}, \tau_m^2 ) \prod_{v=1}^{S_m} \text{N}(\lambda_{u,v}^{(m)} \mid 0, 1) \prod_{v=S_m+1}^{V_m} \text{N}(\lambda_{u,v}^{(m)} \mid \lambda_{u,v-S_m}^{(m)}, \omega_m^2) \right) \end{split} \end{equation}\]
\[\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}\]
RWSeas(n, s = 1, s_seas = 1, along = NULL)n is \(S_m\)s is \(A_{\tau}^{(m)}\)s_seas is \(A_{\omega}^{(m)}\)along used to identify “along” and “by” dimensions\[\begin{align} \beta_{u,v}^{(m)} & = \alpha_{u,v} + \lambda_{u,s_v}, \quad u = 1, \cdots, U_m, \quad v = 1, \cdots, V_m \\ \alpha_{u,v}^{(m)} & \sim \text{N}(0, 1), \quad u = 1, \cdots, U_m, \quad \quad v = 1,2, \\ \alpha_{u,v}^{(m)} & \sim \text{N}(2 \alpha_{u,v-1}^{(m)} - \alpha_{u,v-2}^{(m)}, \tau_m^2), \quad u = 1, \cdots, U_m, \quad v = 3, \cdots, V_m \\ \lambda_{u,v}^{(m)} & \sim \text{N}(0, 1), \quad u = 1, \cdots, U_m, \quad v = 1, \cdots, S_m \\ \lambda_{u,v}^{(m)} & \sim \text{N}(\lambda_{u,v-S_m}^{(m)}, \omega_m^2), \quad u = 1, \cdots, U_m, \quad v = S_m + 1, \cdots, V_m \tag{4.1} \\ \\ \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}\]
We allow \(A_{\omega}^{(m)2}\) to be set to zero, in which case (4.1) reduces to \[\begin{equation} \lambda_{u,v}^{(m)} = \lambda_{u,v-S_m}^{(m)}, \quad u = 1, \cdots, U_m, \quad v = S_m + 1, \cdots, V_m, \end{equation}\] implying that seasonal effects are constant across years.
\[\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, 1 ) \prod_{v=2}^{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} \text{N}(\lambda_{u,v}^{(m)} \mid 0, 1) \prod_{v=S_m+1}^{V_m} \text{N}(\lambda_{u,v}^{(m)} \mid \lambda_{u,v-S_m}^{(m)}, \omega_m^2) \right) \end{split} \end{equation}\]
\[\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}\]
RW2Seas(n, s = 1, s_seas = 1, along = NULL)n is \(S_m\)s is \(A_{\tau}^{(m)}\)s_seas is \(A_{\omega}^{(m)}\)along used to identify “along” and “by” dimensions\[\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 u = 1, \cdots, U_m, \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 individual \(\phi_k^{(m)}\) is restricted to the interval \((-1, 1)\), and the \(\phi_k^{(m)}\) are jointly restricted to values that yield stationary models. Let \[\begin{equation} r^{(m)} = \sqrt{\phi_1^{(m)2} + \cdots + \phi_{K_m}^{(m)2}}. \end{equation}\] We assign \(r^{(m)}\) the prior \[\begin{equation} r^{(m)} \sim \text{Beta}(2, 2). \end{equation}\]
\[\begin{equation} \text{N}^+\left(\tau_m \mid 0, A_{\tau}^{(m)2} \right) \text{Beta}\left( r^{(m)} \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.
\[\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}\]
AR(n = 2, s = 1, along = NULL)n is \(K_m\)s is \(A_{\tau}^{(m)}\)along is used to indentify the “along” and “by” dimensionsSpecial case or AR, but with extra options for autocorrelation coefficient.
\[\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}(2, 2) \\ \tau_m & \sim \text{N}^+\left(0, A_{\tau}^{(m)2}\right). p\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\).
\[\begin{equation} \text{N}(\tau_m \mid 0, A_{\tau}^{(m)2}) \text{Beta}( \phi_m^{\prime} \mid 2, 2) \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}\]
\[\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}\]
AR1(min = 0.8, max = 0.98, s = 1, along = NULL)min is \(a_{0m}\)max is \(a_{1m}\)s is \(A_{\tau}^{(m)}\). Defaults to 1.along is used to identify “along” and “by” dimensionsThe defaults for min and max are based on the defaults for function ets() in R package forecast (Hyndman and Khandakar 2008).
\[\begin{align} \beta_{u,v}^{(m)} & \sim \text{N}(\alpha_u^{(m)} + v \eta_u^{(m)}, \tau_m^2), \quad u = 1,\cdots,U_m, \quad v = 1, \cdots, V_m \\ \alpha_u^{(m)} & \sim \text{N}(0, 1), \quad u = 1,\cdots,U_m \\ \eta_u^{(m)} & \sim \text{N}\left(0, A_{\eta}^{(m)2}\right), \quad u = 1, \cdots, U_m \\ \tau_m & \sim \text{N}^+\left(0, A_{\tau}^{(m)2}\right) \end{align}\]
\[\begin{equation} \text{N}(\tau_m \mid 0, A_{\tau}^{(m)2}) \prod_{u=1}^{U_m} \text{N}(\alpha_u^{(m)} \mid 0, 1) \text{N}(\eta_u^{(m)} \mid 0, A_{\eta}^{(m)2}) \prod_{u=1}^{U_m} \prod_{v=1}^{V_m} \text{N}\left(\beta_{u,v}^{(m)} | \alpha_u^{(m)} + v \eta_u^{(m)}, \tau_m^2 \right) \end{equation}\]
\[\begin{equation} \beta_{u,V_m + h}^{(m)} \sim \text{N}(\alpha_u^{(m)} + (V_m + h) \eta_u^{(m)}, \tau_m^2) \end{equation}\]
Lin(s = 1, sd = 1, along = NULL)s is \(A_{\tau}^{(m)}\)sd is \(A_{\eta}^{(m)}\)along is used to indentify “along” and “by” dimensions\[\begin{align} \beta_{u,v}^{(m)} & = \alpha_u^{(m)} + \eta_u^{(m)} v + \epsilon_{u,v}^{(m)}, \quad u = 1, \cdots, U_m, \quad v = 1, \cdots, V_m \\ \alpha_u^{(m)} & \sim \text{N}\left(0, 1\right), \quad u = 1, \cdots, U_m \\ \eta_u^{(m)} & \sim \text{N}\left(0, A_{\eta}^{(m)2}\right), \quad u = 1, \cdots, U_m \\ \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 u = 1, \cdots, U_m, \quad v = K_m + 1, \cdots, V_m. \end{align}\]
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)}\) is restricted to the interval \((-1, 1)\), and the \(\phi_k^{(m)}\) are jointly restricted to values that yield stationary models. Let \[\begin{equation} r^{(m)} = \sqrt{\phi_1^{(m)2} + \cdots + \phi_{K_m}^{(m)2}}. \end{equation}\] We assign \(r^{(m)}\) the prior \[\begin{equation} r^{(m)} \sim \text{Beta}(2, 2). \end{equation}\]
\[\begin{equation} \text{N}^+\left(\tau_m \mid 0, A_{\tau}^{(m)2} \right) \text{Beta}\left( r_k^{(m)} \mid 2, 2 \right) \prod_{u=1}^{U_m} \text{N}(\alpha_u^{(m)} \mid 0, 1) \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{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.
\[\begin{align} \beta_{u, V_m + h}^{(m)} & = \alpha_u^{(m)} + \eta_u^{(m)} (V_m + h) + \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}\]
Lin_AR(s = 1, sd = 1, along = NULL)s is \(A_{\tau}^{(m)}\)sd is \(A_{\eta}^{(m)}\)along is used to indentify “along” and “by” variables\[\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 4.4).
\(\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\).
\[\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}\]
Terms with a P-Spline prior cannot be forecasted.
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” variablesAge 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 11.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}\]
\[\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
Terms with an SVD prior cannot be forecasted.
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.
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,
\[\begin{align} \pmb{\beta}_{u,t}^{(m)} & = \pmb{F}^{(m)} \pmb{\alpha}_{u,t}^{(m)} + \pmb{g}^{(m)}, \quad u = 1, \cdots, U_m; \quad t = 1, \cdots, T \\ \alpha_{u,k,1}^{(m)} & \sim \text{N}(0, 1), \quad u = 1, \cdots, U_m, k = 1, \cdots, K_m \\ \alpha_{u,k,t}^{(m)} & \sim \text{N}(\alpha_{u,k,t-1}^{(m)}, \tau_m^2), \quad u = 1, \cdots, U_m, k = 1, \cdots, K_m; t = 2, \cdots, T \\ \tau_m & \sim \text{N}^+\left(0, A_{\tau}^{(m)2}\right) \end{align}\]
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, 1) \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, 1) \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}\]
TODO - write
\[\begin{align} \alpha_{u,k,T+h}^{(m)} & \sim \text{N}(\alpha_{u,k,T+h-1}^{(m)}, \tau_m^2), \quad u = 1, \cdots, U_m; \quad k = 1, \cdots, K_m \\ \pmb{\beta}_{u,T+h}^{(m)} & = \pmb{F}^{(m)} \pmb{\alpha}_{u,T+h}^{(m)} + \pmb{g}^{(m)}, \quad u = 1, \cdots, U_m \end{align}\]
SVD_RW(ssvd, n_comp = NULL, s = 1, indep = TRUE)where
- ssvd is an object containing \(\pmb{F}\) and \(\pmb{g}\)
- n_comp is \(K_m\) (which defaults to ceiling(n/2), where n is the number of components in ssvd
- s is \(A_{\tau}^{(m)}\)
- indep determines whether and independent or joint model will be used if the term being modelled contains a sex or gender variable.
Same structure as SVD_RW(). TODO - write details.
Same structure as SVD_RW(). TODO - write details.
Same structure as SVD_RW(). TODO - write details.
Elements of \(\pmb{\beta}^{(m)}\) are treated as known with certainty.
Known priors make no contribution to the posterior density.
Main effects with a known prior cannot be forecasted.
Known(values)values is a vector containing the \(\beta_j^{(m)}\).The columns of matrix \(\pmb{Z}\) are assumed to be standardised to have mean 0 and standard deviation 1. \(\pmb{Z}\) does not contain a column for an intercept.
We implement two priors for coefficient vector \(\pmb{\zeta}\). The first prior is designed for the case where \(P\), the number of colums of \(\pmb{Z}\), is small, and most \(\zeta_p\) are likely to distinguishable from zero. The second prior is designed for the case where \(P\) is large, and only a few \(\zeta_p\) are likely to be distinguishable from zero.
\[\begin{align} \zeta_p \mid \varphi & \sim \text{N}(0, \varphi^2) \\ \varphi & \sim \text{N}^+(0, 1) \end{align}\]
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:
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).
TODO - write
TODO - write
set_covariates(formula, data = NULL, n_coef = NULL)formula is a one-sided R formula describing the covariates to be useddata A data frame. If a value for data is supplied, then formula is interpreted in the context of this data frame. If a value for data is not supplied, then formula is interpreted in the context of the data frame used for the original call to mod_pois(), mod_binom(), or mod_norm().n_coef is the effective number of non-zero coefficients. If a value is supplied, the shrinkage prior is used; otherwise the standard prior is used.Examples:
set_covariates(~ mean_income + distance * employment)
set_covariates(~ ., data = cov_data, n_coef = 5)Use exponential distribution, parameterised using mean, \[\begin{equation} \xi \sim \text{Exp}(\mu_{\xi}) \end{equation}\]
\[\begin{equation} p(\xi) = \frac{1}{\mu_{\xi}} \exp\left(\frac{-\xi}{\mu_{\xi}}\right) \end{equation}\]
set_disp(mean = 1)mean is \(\mu_{\xi}\)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:
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}\]
Running TMB yields a set of means \(\pmb{m}\), and a precision matrix \(\pmb{Q}^{-1}\), which together define the approximate joint posterior distribution of
We use \(\tilde{\pmb{\theta}}\) to denote a vector containing all these quantities.
We perform a Cholesky decomposition of \(\pmb{Q}^{-1}\), to obtain \(\pmb{R}\) such that \[\begin{equation} \pmb{R}^{\top} \pmb{R} = \pmb{Q}^{-1} \end{equation}\] We store \(\pmb{R}\) as part of the model object.
We draw generate values for \(\tilde{\pmb{\theta}}\) by generating a vector of standard normal variates \(\pmb{z}\), back-solving the equation \[\begin{equation} \pmb{R} \pmb{v} = \pmb{z} \end{equation}\] and setting \[\begin{equation} \tilde{\pmb{\theta}} = \pmb{v} + \pmb{m}. \end{equation}\]
Next we convert any transformed hyper-parameters back to the original units, and insert values for \(\pmb{\beta}^{(m)}\) for terms that have Known priors. We denote the resulting vector \(\pmb{\theta}\).
Finally we draw from the distribution of \(\pmb{\gamma} \mid \pmb{y}, \pmb{\theta}\) using the methods described in Section 2.
NOTE - We are ignoring covariates for the moment. We can probably just subtract the covariates term from \(\pmb{\mu}\) and then proceed as before.
Although the sum \(\pmb{\mu} = \sum_{m=0}^{M} \pmb{X}^{(m)} \pmb{\beta}^{(m)}\) is well identified by the data, the individual \(\pmb{\beta}^{(m)}\) are not. For instance, adding \(\Delta\) to \(\pmb{\beta}^{(0)}\) and subtracting it from every term in \(\pmb{\beta}^{(m)}\) for some \(m > 0\) leaves the likelihood unchanged. The use of proper priors for the \(\pmb{\beta}^{(m)}\) mean that posterior distributions for the \(\pmb{\beta}^{(m)}\), and for related terms such as trend, cyclical, and seasonal effects, are proper. However, they can be diffuse and hard to interpret.
To help with interpretation of the \(\pmb{\beta}^{(m)}\) and related terms, we standardize. We in fact appy two forms of standardization, each designed for a different purpose.
The first type of standardization is applied only to the \(\pmb{\beta}^{(m)}\). The aim is to partition the total variation in \(\mu_i\) into variation associated with each main effect or interaction, in the style, for instance, of Section 15.6 of Gelman et al. (2014). Let \(\tilde{\pmb{\beta}}^{(m)}\) be the standardized version of \(\pmb{\beta}^{(m)}\). We obtain the \(\tilde{\pmb{\beta}}^{(m)}\) through the following algorithm:
Set \[\begin{equation} \pmb{r}^{(0)} = \pmb{\mu} \end{equation}\]
For \(m = 0, \cdots, M\):
Term-level standardization aims to clarify the behaviour of the terms making up the model, including trend, cyclical, and seasonal effects. The intercept term is left untouched. All other \(\pmb{\beta}^{(m)}\) are centered across the “along” dimension, and each of the “by” dimensions.
The forecasted values for the time-varying \(\pmb{\beta}^{(m)}\) are shifted up or down so that they line up with the estimated values. The standardization algorithm is then applied to these shifted values.
The forecasted values for the time-varying \(\pmb{\beta}^{(m)}\), and for SVD coefficients, and trend, cyclical, and seasonal components, are shifted up or down so that they line up with the estimated values. All terms are then centered along all dimensions other than time.
To generate one set of simulated values, we start with values for exposure, trials, or weights, \(\pmb{w}\), and possibly covariates $, then go through the following steps:
\[\begin{equation} y_i^{\text{rep}} \sim \text{Poisson}(\gamma_i w_i) \end{equation}\]
\[\begin{align} y_i^{\text{rep}} & \sim \text{Poisson}(\gamma_i^{\text{rep}} w_i) \\ \gamma_i^{\text{rep}} & \sim \text{Gamma}(\xi^{-1}, (\xi \mu_i)^{-1}) \end{align}\] which is equivalent to \[\begin{equation} y_i^{\text{rep}} \sim \text{NegBinom}\left(\xi^{-1}, (1 + \mu_i w_i \xi)^{-1}\right) \end{equation}\]
\[\begin{equation} y_i^{\text{rep}} \sim \text{Binomial}(w_i, \gamma_i) \end{equation}\]
\[\begin{align} y_i^{\text{rep}} & \sim \text{Binomial}(w_im \gamma_i^{\text{rep}}) \\ \gamma_i^{\text{rep}} & \sim \text{Beta}\left(\xi^{-1} \mu_i, \xi^{-1}(1 - \mu_i)\right) \end{align}\]
\[\begin{equation} y_i^{\text{rep}} \sim \text{N}(\gamma_i, \xi^2 / w_i) \end{equation}\]
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.
replicate_data(x, condition_on = c("fitted", "expected"), n = 20)| Quantity | Definition |
|---|---|
| \(i\) | Index for cell, \(i = 1, \cdots, n\). |
| \(y_i\) | Value for outcome variable. |
| \(w_i\) | Exposure, number of trials, or weight. |
| \(\gamma_i\) | Super-population rate, probability, or mean. |
| \(\mu_i\) | Cell-specific mean. |
| \(\xi\) | Dispersion parameter. |
| \(g()\) | Log, logit, or identity function. |
| \(m\) | Index for intercept, main effect, or interaction. \(m = 0, \cdots, M\). |
| \(j\) | Index for element of a main effect or interaction. |
| \(u\) | Index for combination of ‘by’ variables for an interaction. \(u = 1, \cdots U_m\). \(U_m V_m = J_m\) |
| \(v\) | Index for the ‘along’ dimension of an interaction. \(v = 1, \cdots V_m\). \(U_m V_m = J_m\) |
| \(\beta^{(0)}\) | Intercept. |
| \(\pmb{\beta}^{(m)}\) | Main effect or interaction. \(m = 1, \cdots, M\). |
| \(\beta_j^{(m)}\) | \(j\)th element of \(\pmb{\beta}^{(m)}\). \(j = 1, \cdots, J_m\). |
| \(\pmb{X}^{(m)}\) | Matrix mapping \(\pmb{\beta}^{(m)}\) to \(\pmb{y}\). |
| \(\pmb{Z}\) | Matrix of covariates. |
| \(\pmb{\zeta}\) | Parameter vector for covariates \(\pmb{Z}^{(m)}\). |
| \(A_0\) | Scale parameter in prior for intercept \(\beta^{(0)}\). |
| \(\tau_m\) | Standard deviation parameter for main effect or interaction. |
| \(A_{\tau}^{(m)}\) | Scale parameter in prior for \(\tau_m\). |
| \(\pmb{\alpha}^{(m)}\) | Parameter vector for P-spline and SVD priors. |
| \(\alpha_k^{(m)}\) | \(k\)th element of \(\pmb{\alpha}^{(m)}\). \(k = 1, \cdots, K_m\). |
| \(\pmb{V}^{(m)}\) | Covariance matrix for multivariate normal prior. |
| \(h_j^{(m)}\) | Linear covariate |
| \(\eta^{(m)}\) | Parameter specific to main effect or interaction \(\pmb{\beta}^{(m)}\). |
| \(\eta_u^{(m)}\) | Parameter specific to \(u\)th combination of ‘by’ variables in interaction \(\pmb{\beta}^{(m)}\). |
| \(A_{\eta}^{(m)}\) | Standard deviation in normal prior for \(\eta_m\). |
| \(\omega_m\) | Standard deviation of parameter \(\eta_c\) in multivariate priors. |
| \(\phi_m\) | Correlation coefficient in AR1 densities. |
| \(a_{0m}\), \(a_{1m}\) | Minimum and maximum values for \(\phi_m\). |
| \(\pmb{B}^{(m)}\) | B-spline matrix in P-spline prior. |
| \(\pmb{b}_k^{(m)}\) | B-spline. \(k = 1, \cdots, K_m\). |
| \(\pmb{F}^{(m)}\) | Matrix in SVD prior. |
| \(\pmb{g}^{(m)}\) | Offset in SVD prior. |
| \(\pmb{\beta}_{\text{trend}}\) | Trend effect. |
| \(\pmb{\beta}_{\text{cyc}}\) | Cyclical effect. |
| \(\pmb{\beta}_{\text{seas}}\) | Seasonal effect. |
| \(\varphi\) | Global shrinkage parameter in shrinkage prior. |
| \(A_{\varphi}\) | Scale term in prior for \(\varphi\). |
| \(\vartheta_p\) | Local shrinkage parameter in shrinkage prior. |
| \(p_0\) | Expected number of non-zero coefficients in \(\pmb{\zeta}\). |
| \(\hat{\sigma}\) | Empirical scale estimate in prior for \(\varphi\). |
| \(\pi\) | Vector of hyper-parameters |
Let \(\pmb{A}\) be a matrix of age-specific estimates from an international database, transformed to take values in the range \((-\infty, \infty)\). Each column of \(\pmb{A}\) represents one set of age-specific estimates, such as log mortality rates in Japan in 2010, or logit labour participation rates in Germany in 1980.
Let \(\pmb{U}\), \(\pmb{D}\), \(\pmb{V}\) be the matrices from a singular value decomposition of \(\pmb{A}\), where we have retained the first \(K\) components. Then \[\begin{equation} \pmb{A} \approx \pmb{U} \pmb{D} \pmb{V}. \tag{11.1} \end{equation}\]
Let \(m_k\) and \(s_k\) be the mean and sample standard deviation of the elements of the \(k\)th row of \(\pmb{V}\), with \(\pmb{m} = (m_1, \cdots, m_k)^{\top}\) and \(\pmb{s} = (s_1, \cdots, s_k)^{\top}\). Then \[\begin{equation} \tilde{\pmb{V}} = (\text{diag}(\pmb{s}))^{-1} (\pmb{V} - \pmb{m} \pmb{1}^{\top}) \end{equation}\] is a standardized version of \(\pmb{V}\).
We can rewrite (11.1) as \[\begin{align} \pmb{A} & \approx \pmb{U} \pmb{D} (\text{diag}(\pmb{s}) \tilde{\pmb{V}} + \pmb{m} \pmb{1}^{\top}) \\ & = \pmb{F} \tilde{\pmb{V}} + \pmb{g} \pmb{1}^{\top}, \tag{11.2} \end{align}\] where \(\pmb{F} = \pmb{U} \pmb{D} \text{diag}(\pmb{s})\) and \(\pmb{g} = \pmb{U} \pmb{D} \pmb{m}\).
Let \(\tilde{\pmb{v}}_l\) be a randomly-selected column from \(\tilde{\pmb{V}}\). From the construction of \(\tilde{\pmb{V}}\), and the orthogonality of the rows of \(\pmb{V}\), we have \(\text{E}[\tilde{\pmb{v}}_l] = \pmb{0}\) and \(\text{var}[\tilde{\pmb{v}}_l] = \pmb{I}\). This implies that if \(\pmb{z}\) is a vector of standard normal variables, then \[\begin{equation} \pmb{F} \pmb{z} + \pmb{g} \end{equation}\] should look approximately like a randomly-selected column from the original data matrix \(\pmb{A}\).