Given the sparse nature of the available data for this SAE analysis, we used the Bayesian two-stage logistic normal (TSLN) approach we proposed recently [21]. Our previous study showed that the TSLN approach could outperform commonly used area [17, 18, 54] and individual level [55] models both in a simulation study focusing on sparse survey data and an application using the 2017–18 NHS data. The two-stage structure of the TSLN approach includes an individual level stage 1 model, followed by an area level stage 2 model.

The same TSLN approach, with very similar components, was chosen to be applied to all eight risk factor measures. The selection of fixed and random effect structures for the two models was guided by the goal of achieving a balance between parsimony across risk factor measures and predictive performance. We followed the advice by Goldstein [56] and initially used frequentist algorithms to select fixed and random effects, with fully Bayesian inference via Markov chain Monte Carlo (MCMC) for final model checking. Further details regarding model selection are given in Section E of the Additional File 1.

Let \(y_ \in \\) be the binary value from the NHS for sampled individual \(j = 1, \dots , n_i\) in SA2 \(i = 1, \dots , m\), where \(n_i\) is the sample size in SA2 i. Further, let \(m = 1694\) and \(M = 2221\) be the number of sampled and total number of SA2s, respectively. The goal of this analysis is to generate estimates of the true proportions of each risk factor measure, \(\varvec = \left( \mu _1, \dots , \mu _M \right)\).

In this analysis, we used two versions of the survey weights, \(w^}_\), provided by the ABS [55, 57] to correct for sampling bias and promote design-consistency. The first, \(w_\), was used for direct estimation and the second, \(\tilde_\), was used in the stage 1 model (see Section C.1 in the Additional File 1). Using the survey weights, small area proportion estimates can be computed using the Hajek [58] direct estimator,

$$\begin \hat^D_i = \frac^ w_ y_ }, \end$$

(1)

with an approximate sampling variance of [54, 59],

$$\begin & \psi _i^D = \widehat} \left( \hat_i^D \right) = \frac \left( 1 - \frac \right) \left( \frac \right) \nonumber \\ & \quad \sum _^ \left( w_^2 \left( y_ - \hat_i^D \right) ^2 \right) . \end$$

(2)

Direct estimators, such as Eqs. (1) and (2), have low variance and are design-unbiased for \(\mu _i\) when \(n_i\) is large, but have high variance when \(n_i\) is small [13].

The stage 1 model is a Bayesian pseudo-likelihood logistic mixed model. Let \(\pi _\) be the probability of \(y_ = 1\) for sampled individual j in SA2 i. Following the notation of Parker et al. [55], we represent the pseudo-likelihood for a probability density, p(.), as \(p\left( y_ \right) ^_}\). Pseudo-likelihood is used to ensure the predictions from the logistic model are approximately unbiased under the sample design [60, 61]. Thus, the stage 1 model likelihood is given by,

$$\begin y_ \sim \text \left( \pi _ \right) ^_}, \end$$

(3)

where \(\text \left( \pi _ \right)\) is modelled using a generic linear predictor that is application-specific. In this work, we used several unique components summarised in Fig. 2. The linear predictor included eight individual level categorical covariates and seven area level covariates as fixed effects. Unstructured individual and SA2 level random effects were also applied. In addition to these, borrowing ideas from MrP [62], we included two hierarchical random effects based on categorical covariates that were themselves derived from the interaction of numerous individual level demographic and health covariates. A discussion of the priors used is given on the subsequent page. More details can be found in Section C of the Additional File 1.

Schematic describing the components of the linear predictor for \(\text \left( \pi _ \right)\) in the stage 1 model. *The non-outcome risk factor categorical covariate was derived from the interaction of the binary risk factor outcomes not directly associated with the risk factor being modelled. For more details see Section C of the Additional File 1. SA2 Statistical area level 2

After fitting the stage 1 model, the individual level predictions are aggregated to the area level, producing stage 1 (S1) proportion estimates \(\hat^, (t)}_i\) using Eq. (1), and sampling variances, \(\psi ^, (t)}_i = \widehat} \left( \hat_i^D \right) + \widehat} \left( \hat^_i \right)\), for all posterior MCMC draws, \(t=1, \dots , T\) [21], where the function to compute the sampling variance, \(\widehat}(.)\), is given in Eq. (2) and \(\hat^_i = n_i^ \left( \sum _^ w_ \left( \pi ^_ - y_ \right) \right)\) quantifies the level of smoothing achieved by using \(\pi _\) instead of \(y_\).

Using the common logistic transformation [18, 54], let

$$\begin \hat_i^, (t)}= & \text \left( \hat_i^, (t)} \right) \end$$

(4)

$$\begin \tau _i^, (t)}= & \psi _i^, (t)} \left[ \hat_i^, (t)} \left( 1 - \hat_i^, (t)} \right) \right] ^, \end$$

(5)

thereby permitting the use of a Gaussian likelihood in the second stage model. Let \(\bar_i^}\) be the empirical posterior mean of \(\tau _i^}\) and \(\widehat} \left( \hat_i^} \right)\) be the empirical posterior variance of \(\hat_i^}\). Finally, by selecting a random subset of the posterior draws, say \(\widetilde\), let \(\hat}^}_i = \left( \hat^, (1)}_i, \dots , \hat^, (\widetilde)}_i \right)\).

The stage 2 model is a Bayesian spatial Fay-Herriot [14] model. Unlike previous two-stage approaches [26, 27], we accommodate some of the uncertainty inherent in fitting the stage 1 model by using the vector \(\hat}^}_i\) as input to the stage 2 model. The stage 2 model likelihood for the posterior draws from the stage 1 model is,

$$\begin \hat}_i^} \sim \text \left( \theta _i , \bar_i^} + \widehat} \left( \hat_i^} \right) \right) \end$$

(6)

where \(\theta _i\) is modelled using a generic linear predictor that is problem specific. The final proportion/prevalence estimate for the ith SA2, denoted \(\mu _i\), is given by the posterior distribution of \(\text ^ \left( \theta _i \right)\). To ensure that posterior uncertainty remains unaffected by the choice of \(\widetilde\), we downscale the likelihood contribution by \(1/\widetilde\).

In this work, we used several unique components for the linear predictor of \(\theta _i\) which are summarised in Fig. 3. The linear predictor included the SES index deciles and remoteness as standard fixed effects. In addition, PC1 to PC6 were used as fixed effects with coefficients varying according to remoteness. The linear predictor also included an external latent field constructed from the SHAA’s estimates and a BYM2 spatial random effect [63] at the SA2 level. Given we did not include SA3 level census covariates, an unstructured random effect at the SA3 level was employed. To smooth unstable variances we used the generalized variance function [12, 64, 65] described in Section C.4.6 of the Additional File 1. More details can be found in Section C of the Additional File 1.

Schematic describing the components of the linear predictor for \(\theta _i\) in the stage 2 model. For more details see Section C of the Additional File 1. SA2: Statistical area level 2; SA3: Statistical area level 3; SHAA: Social Health Atlases of Australia

The Bayesian models described above are completed by the specification of priors. Given the complexity of the two models, in this work generic weakly informative priors were adopted based on preliminary analysis of the data [66]. In both models, all fixed effect coefficients were given \(\text \left( 0, 2^2 \right)\) priors with intercepts given a student-\(t\left( 0, 2^2, \text = 3 \right)\). We used \(\text ^\left( 0, 1^2 \right)\) and \(\text ^\left( 0, 2^2 \right)\) priors for all standard deviation terms in the stage 1 and stage 2 models, respectively. The mixing parameter in the BYM2 [63] random effect was given a \(\text \left( 0,1 \right)\) prior (see Section C of the Additional File 1).

We conducted sensitivity analysis by using more, \(\text \left( 0, 1^2 \right)\), and less, \(\text \left( 0, 100^2 \right)\), informative priors for the fixed effects in both models. We also experimented with using exponential priors with rates of 0.5 and 1 for standard deviation terms. Finally, we examined model fit when using an informative Beta prior for the mixing parameter. We found that the model fit and prevalence estimates were unaffected by these prior changes. The chosen priors gave superior sampling efficiency and convergence.

For validation of the small area estimates, we adopted a dual approach, using both internal and external methods. See Section C.5 in the Additional File 1 for details.

Internal validation involved a fully Bayesian benchmarking procedure [29] that adjusts the results obtained in the stage 2 model by penalizing discrepancies between modelled and direct estimates. Unlike previous benchmarking approaches that adjust the point estimates only [13, 67], Bayesian benchmarking adjusts the entire posterior — automatically accounting for benchmarking-induced uncertainty.

In this work we simultaneously enforced two benchmarks referred to as “state” and “major-by-state”. The state benchmark had seven groups which were composed of the states and territories of Australia (except the Northern Territory, which was not benchmarked due to ABS instruction [57]).

The major-by-state benchmark had 12 groups, composed of the interaction of the states and territories of Australia (except the Northern Territory) and dichotomous remoteness (major city vs non-major city). Thus, for each state, apart from Tasmania (where all areas were non-major city), and the Australian Capital Territory (where 96% of areas were major city), each SA2 was benchmarked differently depending on whether the area was in a major city or not.

External validation was performed by comparing the estimates to those from the SHAA at the PHA level and the overall trends observed in the modelled results with the general findings from other Australian health surveys conducted on specific sub-populations, such as states [68] or First Nations Australians [69]. Although this validation affirmed the validity and reliability of our estimates in general, it was particularly helpful in assessing the credibility of estimates for areas that could not be benchmarked.

We used fully Bayesian inference using MCMC via the R package rstan Version 2.26.11 [70]. Where possible we used the non-mean centered parameterization for random effects and the QR decomposition for fixed effects [71]. The stan code for the stage 1 and stage 2 models can be found on GitHub [72].

For the stage 1 model we used 1000 warmup and 1000 post-warmup draws for each of the four chains, feeding a random subset of 500 posterior draws from the stage 1 to the stage 2 model. For the stage 2 model we used 3000 warmup and 3000 post-warmup draws for each of four chains. For storage reasons we thinned the final posterior draws from the stage 2 model by 2, resulting in 6000 useable posterior draws.

Convergence of the models was assessed using trace and autocorrelation plots, effective sample size and \(\hat\) [73]. While convergence ranged slightly between risk factors, all the proportion parameters, \(\varvec = \left( \mu _1, \dots , \mu _M \right)\), had \(\hat < 1.03\), with 96% having effective sample sizes \(>1000\) and 99% having \(\hat < 1.01\).

Estimates from the benchmarked stage 2 model were reported in a variety of forms, including absolute, relative and classification measures. For point estimates we used posterior medians and for uncertainty intervals we used 95% highest posterior density intervals (HPDIs). We used the modelled proportions as the absolute indicator and odds ratios (ORs) as the relative indicator. The ORs for the tth posterior draw were derived as,

$$\begin \text ^_i=\, & \frac_i/(1-\mu ^_i)}^D/(1-\hat^D)} \end$$

(7)

with \(\hat^D\) being the national prevalence estimate for the risk factor measure. An OR above one indicates that the SA2 has a prevalence higher than the national average.

In addition to using point estimates and credible intervals to summarize the ORs, we also used the exceedence probability (EP) [31, 53, 74].

$$\begin EP_i = \frac \sum _t \mathbb \left( \text ^_i > 1 \right) \end$$

(8)

Generally an EP above 0.8 (or below 0.2) is considered to provide evidence that the proportion in the corresponding SA2 was substantially higher (or lower) than the national average, respectively [75]. Note that the exceedance probabilities calculated using either ORs or prevalence are identical.

To facilitate decision-making, we classified SA2s by assessing whether their individual and neighbor values (i.e. clusters [76, 77]) were significantly different to the national average. In this work, these classifications were called evidence classifications. Any area classified as HC, H, L, or LC has an exceedance probability suggesting that the modelled prevalence is significantly different to the national average; HC or H denotes higher, while L or LC denotes lower. The difference between HC and H (or LC and L) is that the former provides an indication of clustering of areas, while the later only indicates significance of the area itself. If an area is not classified according to the criteria above (defined as None (“N”)) the modelled estimate is not sufficiently different to the national average. See details in Section D.3 of the Additional File 1.

Code to produce subsequent plots is available on GitHub [72].