Bayesian Sample Size Determination for Planning Hierarchical Bayes Small Area Estimates
Peter Dutey-Magni

	A superconvergent HDG method for Stokes flow with strongly enforced symmetry of the stress tensor
    
			  Matteo Giacomini

             Laboratori de Càlcul Numèric (LaCàN), 

             ETS de Ingenieros de Caminos, Canales y Puertos, 

             Universitat Politècnica de Catalunya, Barcelona, Spain

             Alexandros Karkoulias

			  Laboratori de Càlcul Numèric (LaCàN), 

             ETS de Ingenieros de Caminos, Canales y Puertos, 

             Universitat Politècnica de Catalunya, Barcelona, Spain

             Centre Internacional de Metodes Numerics en Enyinyeria (CIMNE), 

             Campus Nord UPC, Barcelona, Spain

			  Ruben Sevilla

             Zienkiewicz Centre for Computational Engineering, 

             College of Engineering, Swansea University, Wales, UK 

             Antonio Huerta

             Laboratori de Càlcul Numèric (LaCàN), 

             ETS de Ingenieros de Caminos, Canales y Puertos, 

             Universitat Politècnica de Catalunya, Barcelona, Spain


    December 30, 2023
============================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================




This paper devises a fully Bayesian sample size determination method for hierarchical model-based small area estimation with a decision risk approach. A new loss function specified around a desired maximum posterior variance target implements conventional official statistics criteria of estimator reliability (coefficient of variation of up to 20 per cent). This approach comes with an efficient binary search algorithm identifying the minimum effective sample size needed to produce small area estimates under this threshold constraint. Traditional survey sampling design tools can then be used to plan appropriate data collection using the resulting effective sample size target. This approach is illustrated in a case study on small area prevalence of life limiting health problems for 6 age groups across 1,956 small areas in Northern England, using the recently developed Integrated Nested Laplace Approximation method for spatial generalised linear mixed hierarchical models. 

Keywords: Disease prevalence models, hierarchical Bayes, Integrated Nested Laplace Approximation, official statistics, sample size determination, small area estimation.





§ INTRODUCTION

When small area estimates from survey samples are imprecise, model-based estimation is often used to borrow strength and produce more efficient estimators. Hierarchical Bayes (HB) prediction, in particular, has been shown to yield good results and make efficient use of available data <cit.>. When it comes to routine applications, an question commonly faced by practitioners is how to determine the minimum sample size necessary to achieve a desired precision, generally expressed as a target standard error or RSE applicable to all target small areas. Although some analytical solutions have been proposed in the model-assisted literature, this problem has remained relatively unexamined in model-dependent SAE.

This paper provides a fully Bayesian treatment of this problem. Section <ref> briefly reviews existing work on Bayesian SSD for clinical trials and other investigations, and draw implications for planning model-based SAE studies. Section <ref> describes a model-based simulation procedure to determine ESS requirements under a maximum relative posterior variance constraint. This approach is applied in section <ref> using 2011 UK census data on chronic health and a spatial hierarchical generalised linear mixed model. A design-based simulation study confirms both (a) the validity of solution produced by the SSD algorithm and (b) efficiency savings achieved over traditional survey sampling direct estimation in section <ref>. Finally, section <ref> discusses generalisation to other model-based methods as well as more complex sampling designs.



§ REVIEW OF SSD FOR HIERARCHICAL MODELS

In the frequentist approach, target sample sizes are usually determined by reference to the sampling distribution of the target parameter under a given survey sampling plan. With complex model-assisted or model-dependent statistical designs, this sampling distribution is typically unknown. This is particularly the case when a study examines a model parameter estimated under a multivariate statistical model with covariate inputs.

Amongst those complex statistical designs are two-stage hierarchical models commonly used to estimate some kind of effect in multicentre randomised controlled trials or in multilevel observational studies (e.g. studies of educational attainment across classes and schools). This effect can be a treatment effect, a regression coefficient or some other cluster-specific characteristic. We distinguish between two categories of studies based on their inferential motivation:

	
  * studies aiming to detect a global effect with a predefined statistical power;
	
  * studies aiming to estimate effect sizes within each centre/cluster with a predefined precision.




 §.§ SSD for detection

In category (1), success is defined as collecting a sufficient number of observations across a sufficient number of centres or clusters so that the target statistical power is achieved in the model of interest. Achieving the target statistical power allows study investigators to reject a null hypothesis if the size of the treatment effect exceeds a clinically meaningful threshold set at the design stage. In observational study designs, examples are varied which optimise, for instance, the number of clusters and the number of units within clusters. Several software applications now exist across the educational, behavioural and wider multilevel literature <cit.>. In interventional research, including clinical trials, some design problems have closed-form solutions <cit.>, while others are more complex and require computer-intensive Monte Carlo experiments to reach a solution. <cit.> proposed algorithms relying on a binary sample size search between 0 and the frequentist binomial sample size requirement to optimise the SSD solution at a reasonable computational cost. The adequate sample size is determined as the smallest sample size meeting the chosen constraints.



 §.§ SSD for prediction

The present paper is concerned with studies belonging to category (2), in which success is defined as collecting a sufficient number of observations to estimate the treatment effect with a predefined precision within each centre or cluster. In contrast with category (1), research on SSD for a vector parameter estimated under a working model is not widely developed, with the exception of clinical audit design methodology. In a series of simulation studies, <cit.> explored the effect of the number of clusters and cluster sample sizes on the efficiency of audits under beta-binomial hierarchical models borrowing strength and adjusting for confounders. <cit.> and <cit.> proposed an SSD algorithm in three-stage hierarchical models search for an design for hospital benchmarking clinical audits. In these examples, the chosen SSD criterion is an upper bound for the average width of posterior intervals of centre-level estimates. Their method involves a grid search of the minimum target sample size (a sample size sweep in a predetermined sequence of candidate values). It is implemented with Monte Carlo simulation using the prespecified model, successively: sampling from the model  priors; simulating a realisation from the model; predicting the corresponding parameter; and computing the level of compliance with the target posterior interval width. This type of Monte Carlo solution associated with Bayesian decision theory are attracting growing interest for complex SSD problems in hierarchical designs for medical studies. Another area for application is the design of spatial sampling. Even with the simplest working models, <cit.> conclude that sampling simulations are inevitable, albeit computationally intensive. The authors nevertheless recognise that the latest model estimation methods are likely to make the simulation approach feasible.



 §.§ SSD for small area studies

SSD techniques have scarcely been applied to plan model-based SAE studies, that is when a single working model is used to compute predictions within each small area in a given study population. Although simulation studies are ubiquitous in the SAE literature, they have been used either to illustrate efficiency gains obtained from more complex modelling designs <cit.>, or to validate a small area model against historical data <cit.>.
Closed-form SSD solutions are only available for the simplest models: see <cit.> in the model-assisted literature, and <cit.> who looks at implications of introducing a model covariate on SSD. Just recently, <cit.> proposed solutions for optimal sample allocation across small areas in relation to the Empirical Best Linear Unbiased Prediction. Yet, the availability of SSD methods for real-world statistical needs remains of particular importance with model-based estimation since precision depends not just on the type of predictor selected, but on covariates and random components included in the model <cit.>. The MSE of predictors depends not solely on a known sampling distribution, but also on the working model's covariance structure. To the best of our knowledge, not SSD methodology is available for emerging developments such as spatial smoothing models <cit.>, or small area models borrowing strength from time <cit.>, age or cohort effects <cit.>.



 §.§ Prerequisites of frequentist and Bayesian solutions for SSD

SSD presupposes (1) design constraints in the form of a set of criteria against which to optimise sample design; and (2) prior information on the population of study.

The implementation of design constraints is straightforward in closed-form frequentist solutions which focus on estimator variance targets. In fully Bayesian SSD, design criteria are often treated as a decision rule determined by a loss function <cit.>. A typical decision rule for studies preoccupied with hypothesis testing is based on a function of statistical power (opposite of the false negative rate) and significance level (opposite of the false positive error rate), as illustrated by <cit.>. Yet other rules have been proposed for detection studies; <cit.> considered three such design objectives: achieving a desired average coverage probability for highest posterior density intervals; satisfying a maximum average length for those intervals; and a combination of the two constraints. As for studies interested in prediction, decision rules have in majority been specified in relation to the width of interval estimates; see for instance <cit.> and <cit.>. These overlap conventional government survey design criteria defined in terms of precision, either based on an estimator's margin of error or RSE. For instance the 2000/01 English Local boost of the Labour Force Survey was designed with a frequentist approach to insure an acceptable precision, defined as a maximum RSE of 20 per cent for design-based estimators of economic activity headcounts of 6,000 districts <cit.>. A model-based equivalent is the relative root MSE, while a fully Bayesian equivalent would be to consider a function of the relative posterior variance, for instance the number of  districts which fail to meet a maximum relative posterior variance threshold.

The second prerequisite of SSD, namely prior information around the parameter of interest, is treated very differently in the frequentist and the Bayesian apparach <cit.>. In the frequentist paradigm, SSD is entirely determined by the sampling distribution of the study target parameter, which itself depends on its population variance, which is typically unknown. In the absence of knowledge regarding the population variance of the study parameter, the standard frequentist approach to SSD involves plugging an assumed value for this variance in a closed formula. The outcome is therefore entirely dependent on how conservative this assumption is and it is often necessary to overestimate the population variance.
In contrast, fully Bayesian SSD does not handle unknown parameters using a plug-in method but instead using explicit priors and hyperpriors <cit.>. The literature contains a variety of examples. When determining a sample size for a binomial proportion, <cit.> and <cit.> set the scale and rate (hyperparameters) of the beta distribution (prior) believed to determine an overdispersed binomial distribution of interest. In biomedical research, such hyperpriors are generally elicited from pilot data, previous studies and subjective expert opinion <cit.>.

At the first glance, particularly when relying on pilot or historical data to form a prior, it seems intuitive to use a single prior for both the design and the analysis. In other words, the set of priors used to simulate prediction under various sample size scenarios is also incorporated in the working model. Yet  <cit.> and <cit.>  
have argued in favour of separating design and fitting priors on scientific grounds. Regulations on biomedical trials sometimes impose that the data are analysed under a state of pre-experimental knowledge, that is without incorporating knowledge from data produced in previous studies. Although such historical data can be valuable in optimising SSD, it is not necessarily desirable to introduce them in the analysis itself as this can be left to subsequent meta-analytic studies. Similar constraints may exist in official statistics, where the reliance on informative priors is sometimes subject to objections <cit.>. When the only available source of prior knowledge is historical data, there may be reasons to restrict its use to SSD. This provides further assurance in the solution of SSD while determining sample requirements to produce a sufficiently precise estimate without having to pool data from previous statistical bulletins. This is by no mean the only way to proceed, but it is sometimes desirable to treat the elicitation of design priors and fitting priors separately. 

On the one hand, design priors retain a strong influence on the outcome of any SSD procedure and its success. Particular attention has been given to robust prior elicitation—that is priors that do not convey excessive confidence compared to the existing knowledge, and which are flexible enough to offer protection against misspecified models. 
With regard to SAE, priors can be elicited from previous survey waves or pilots: routine government survey data are typically abundant. In principle, the most simple type of design prior can consist of hyperparameters of random components since they determine the level of shrinkage in HB prediction and, by way of consequence, posterior variance. Such hyperpriors can be elicited from appropriate marginal posteriors obtained from the combination of pilot data with a vague uninformative prior (which can be the fitting prior). Due consideration must be given to how much belief can be placed into the stability of these parameters across years or across surveys, and it may be necessary to apply a small discount to their influence (see ).

On the other hand, fitting priors can be everything between vague and informative. The use of uninformative hyperpriors for random components is not always an option as it often leads to difficulties in estimating models. Both Markov chains Monte Carlo and INLA encounter numerical difficulties with complex models, especially with spatial models.  <cit.> suggest specifying weakly informative priors for Gaussian random effects by using the log Student t distribution (with one or two degrees of freedom) and predefined lower and upper bounds for the range of 95 per cent of realisations. Hyperpriors can then be deduced in the form of the scale and rate of a Gamma distribution. More recently, <cit.> have addressed more complex models with a combination of random effects and proposed a weakly informative `penalised complexity prior' based on some belief of the scale of random effect (standard deviation).





§ MODEL-BASED SSD
 
We consider the estimation of a small area characteristic (e.g. economic status, illness, income, marital status) under a given working model ℳ.
Let the population U be partitioned into small areas d={1,…, D}. 𝐍 = {N_1,…, N_D} and  𝐘 = { Y_1,…, Y_D} respectively denote the population size in area d and the characteristic total in area d: this can be the area headcount of  individuals with the given characteristic, or the area total (such as total income). We are interested in estimating population means Y_d= Y_d/N_d using

	
  * hierarchical working model ℳ, to produce an HB predictor of  Y_d notated θ_d;
	
  * data from auxiliary covariates 𝐗 available for the entire population;
	
  * a sample survey s to be designed.


We seek to determine f, the effective sampling fraction for s, using  Bayesian SSD under some design constraints. We remark that area-specific ESS n_d are such that n_d∼Binomial(N_d, f).



 §.§ Sample size criteria and Bayes decision rule

2.5ptto10pt1.8pt
		q   
			1 j .7 0 0 .7 0 0 cm
			-2 1 m 
			3 2 4 5 4 7  c
			4 9 3 10 2 10 c
			1 10 0 9 0 7  c
			0 3 5 2 5 0 c
			5 -2 2 -3 1 -3 c
			S Q4.5pt

A conventional frequentist criterion of statistical reliability in official statistics is the RSE or CV. A possible Bayesian equivalent is the relative posterior variance of the HB predictor θ_d. Tabular cells of estimates with a relative posterior variance in excess of 20 per cent are to be suppressed from statistical publications. 
It is desirable that the overall rate of cell suppression remains low, for example below a threshold = 0.01. We implement this requirement through a simple design loss function: the proportion of cell suppression in the dataset weighted by the population headcount of those cells. The total loss ℓ(f  |  …) can be thought of as the total headcount of populations eligible to a reliable estimate, whose estimates have to suppressed due to reliability concerns. The weighting introduces a form of trade-off, which prioritises reliable estimates for large populations while being more tolerant of the risk of cell suppression for the smallest cross-classifications, which are inevitably the most demanding in terms of data collection.

    ℓ(f  |  ℳ,  π(τ_γ), π(τ_υ), π(τ_ν), …) =   N^-1∑_d[   N_d I(RSE(θ_d) > 0.2)]


In this expression, I(·) is the indicator function. The overall loss is an intractable function of the sampling fraction f and is conditional on both the working model ℳ and the set of design and fitting priors π(·). Though ℓ has no obvious closed-form expression, it is reasonable to take the premise that it is a monotonically decreasing function of n (for a discussion of the design consistency of HB prediction, see ). This facilitates the evaluation of integral ∫ℓ(f  |  …)   df over a reasonable interval of sampling fractions [a,b].

Many more specifications can be considered for the loss. ℓ can be defined with respect to the width of prediction intervals such as highest posterior density intervals rather than RSE. It can also be restricted to a subset of domains d in the event that quality standards set by statutory or funding requirements do not apply to all domains d. 

Depending on design constraints and costs, it is possible to envisage more sophisticated loss functions inspired from optimal sample design or the `value of information' approach. In particular, it is conceivable to attribute a price to cell suppressions or to penalised ℓ by a marginal cost function of increasing the sampling fraction. These constitute options for more holistic decision rule and lead to an optimum between cost saving and exhaustive publication. Provided that such refinements strictly depend on f and posterior means or variances, they add no further dimension to the SSD equation, and the approach described in this paper should remain entirely applicable.







 §.§ Sample size minimisation

Simulating the relative posterior variance for all possible sample sizes until the minimum quality standard is attained following a sweep approach is computationally cumbersome. We instead opt for a binary search algorithm (see Figures <ref> and <ref>) to minimise ℓ(f). By progressing iteratively towards the solution with steps of decreasing sizes, we considerably reduce the number of simulations. The algorithm starts with f_1, the midpoint of the interval [a,b], and evaluates ℓ(f_1) over many replications of the sampling and model estimation process. This process is then replicated on interval [a, f_1] if ℓ(f_1) can be trusted to fall below the maximum tolerated suppression setting  or [f_1,b] in the opposite case—that is, if more than  per cent of sampling simulations return a strictly positive value for ℓ(f_1). The interval of possible solutions f_k is therefore halved repeatedly until the interval width is small enough to be used; this is determined by a tolerance setting h. The upper bound of the final interval is the target ESS solution. The algorithm can be restarted on the final interval to reach an even narrower h. By construct, the number of iterations required to complete the algorithm  is the smallest integer  k_max such that k_max > (2h)^-1(b-a).

It is worth emphasising that this procedure relies on ℓ being a monotonous (not strictly) decreasing function of sampling fraction f. The algorithm progression is mainly determined by a risk setting γ, the probability that cell suppression will exceed  for any given f_k. The number of simulations L ordered for each sampling fraction f_k should be set in relation to γ; we suggest that the order of magnitude of L should always be at least equal to the order of magnitude of γ^-1.




§ CASE STUDY
 


 §.§ Context

The approach described in section <ref> is illustrated on a typical tabular Office for National Statistics census output. This case study examines sample size requirements and efficiency gains under small area models predicting the age-specific prevalence of  LLTI across small areas in Northern England. The UK government having announced its goal to discontinue full population enumeration through traditional decennial censuses <cit.>, it is expected that small area population health characteristics will in the future have to be estimated using sample surveys and administrative sources. The aim of the present case study is twofold: (i) to apply the SSD procedure set out in section <ref> to reproduce 2011 census health outputs under conventional UK standards of statistical reliability; (ii) to assess efficiency gained from incorporating more advanced hierarchical prior structures into small area models (borrowing strength from neighbouring areas and age groups).


Future census quality standards published by the <cit.> demand a maximum RSE of 20 per cent for population headcounts representing at least 3 per cent of an MSOA's population. Assuming an average resident population of 8,000 across MSOA, this criterion is equivalent to a margin of error of up to ± 94 for a headcount of 240 persons. MSOA zones were originally designed to ensure lower and upper population limits of 5,700 to 11,100 usual residents. By 2014, 10 per cent of these zones no longer met such limits <cit.>. As a consequence, the future Office for National Statistics' census quality standards in effect require a much smaller margin of error for the smallest MSOA.



 §.§ Design constraints on a tabular output

The case study population U is partitioned into MSOA d={1,…, D} and age groups j ={1,…, J}. 𝐍 = {N_11,…, N_JD} and  𝐘 = { Y_11,…, Y_JD} respectively denote headcounts of members of private households of age group i residing in area d and the corresponding number of individuals reporting an LLTI. We are interested in the LLTI prevalence proportions p_jd = Y_jd/N_jd. Table <ref> illustrates the typical tabular census output corresponding to this data. 




In the present case study, the Office for National Statistics' census quality threshold is applied to  tabular cells  Y_jd (headcounts of individuals reporting an LLTI) but not to cells N_jd - Y_jd (headcounts of individuals reporting no LLTI). Although the Office for National Statistics' quality standard would normally apply to these cells as well, we consider that most public health statistics producers would only publish prevalence or disease headcounts estimates, not estimates of disease-free headcounts. For simplicity, we neither predict nor test the precision of LLTI-free cells. Furthermore, we set the tolerable cell suppression rate to = 0. Because we aim for no suppression at all, there are no longer any benefits in weighting the loss function (<ref>) by the underlying population sizes, as discussed in section <ref>. In those conditions, we substitute the loss function used at line <ref> of Algorithm <ref> with the simpler crude proportion of suppression below:

    ℓ(f_k)_l ←∑_j=1^J∑_d=1^D(      I[RSE_jdl > 0.2] )


Since the quality threshold does not apply to all cross-classifications, the above loss ℓ is only computed over cross-classifications jd such that Y_jd/N_d≥0.03. Because Y_jd is in reality unknown, we compute the estimated loss ℓ̂ over cross-classifications jd such that θ_jd/N_d≥0.03. We nevertheless report both the true and the estimated loss in results in order to assess any potential divergence in the algorithm's progression.

Using data from the 2011 UK census table LC3101EWLS <cit.>, 18 age groups were collapsed in order to bring Y_jd/N_d to an average level close to 3 per cent or more. With D=1,956 and setting J=6, we have a total JD=11,736 cells, 49 per cent of which are eligible for the quality standard of a 20 per cent RSE. Other settings are configured at f_a = 0.01, f_b=0.04, h=0.01 and L=100. The acceptable risk of exceeding the tolerable rate of suppression  is set to γ = 0.01.



 §.§ Model specification
  
We consider the below working model (<ref>). Sample counts y_jd are treated as the realisation of a binomial distribution with sample size n_jd and success parameter p_jd=logit^-1(θ_jd). 

    y_jd∼Binomial(logit^-1(θ_jd) , n_jd)      
       θ_jd = β^(1)_j +  X_dβ^(2)_j +  υ_d  + ν_jd

where β^(1) is an i-dimensional vector of age contrasts (with β^(1)_1 = 0 for identifiability); 𝐗 a matrix of area-level scaled covariates; β^(2) is a matrix of coefficients controlling the effect of one standard deviation in the covariates on area-level log-odds of LLTI; a structured ICAR area random effect υ_d; an unstructured (exchangeable) area by age random effect ν_jd.
 
Covariates in 𝐗 are taken from public sources; namely: the MSOA- and district-level ISAR <cit.>;  the 2015 Income Deprivation Affecting Children Index (proportion of all children aged 0–15 years living in income-deprived families based on tax and benefit departmental database, ); the MSOA mean price of residential property sales in 2011 <cit.>; 2011 mortality ratios indirectly standardised by sex and age <cit.>; and the 2011 MSOA Rural Urban Classification consisting of five contrasts <cit.>.
No age-specific area covariates were available, which reduces the predictive capabilities of the model: although ISAR exhibit strong levels of correlation with crude LLTI prevalence across MSOA, associations with age-specific prevalence are weaker.

 

The structured ICAR area effect υ_d is included to model the shared spatial surface in disease prevalence as first introduced by <cit.> under a sum-to-zero constraint for identifiability: ∑_d=1^Dυ_d = 0. Every structured area effect  υ_d is dependent on other area structured effects υ_[d] = {υ_j, j ≠ d }  under the following conditional distribution: 

    υ_d  | υ_[d]∼Normal(∑_d ≠ jw_djυ_jw_dj , 1τ_υ∑_[d] w_dj)


where spatial weights w are taken from the spatial dependence matrix 𝐑_υ^-1 defined as a D × D-dimensional contiguity matrix (Queen's method) as below:

    w_ij  =   
    1          i = j    
    
    1         i, j  are neighbours
    
    0         otherwise


corresponding with the below joint prior density:

    π(υ | τ_υ)   ∝ exp(-τ_υ2 υ' 𝐑_υυ)


ν has a basic normal exchangeable prior centred around zero with a unique precision τ_ν. The sum of υ and ν forms the widely used convolution prior <cit.>. A discrete mixture of normal exchangeable effects with unequal variances was considered to take into account evidence of heteroskedasticity across age groups. This has not led to substantial improvement and has thus been abandoned.

Models are estimated using INLA <cit.>, given the substantial computational gains over Markov Chains Monte Carlo. Simulation results for generalised linear mixed models <cit.> have shown that posterior distributions obtained are virtually perfectly aligned to posteriors obtained with Markov Chains Monte Carlo sampling. Estimation is implemented using software packages  v 1.698,  3.2.1, and package  <cit.> on high performance computer IRIDIS 4, using 16-core nodes, each equipped with 2.6 GHz CPUs and 4 GB RAM.




 §.§ Study and fitting priors

The fitting prior can be formed based on the accumulation of modelling experience on the outcome of interest, particularly around the scale of random effects. It can be difficult to anticipate the magnitude of residual between-area heterogeneity once covariates are used, but the accumulation of evidence over  time can guide priors around the magnitude of area random effects expressed as log-odds. We specify inverse Gamma priors for these effects on the basis of previous evidence that once covariates are introduced, their variance is likely to be near 0.1,
and that half of this between-area heterogeneity is likely to be spatially structured. Following this, we place equal weight on structured and unstructured effect. To achieve a sum variance of 0.1, we need a mean hyperprior precision of 5 for each area effect. Based on plots we also decide that the variance of this hyperprior should be ten times the mean to avoid dominating the data. This leads to hyperpriors Γ(0.5, 0.1) being placed on each of the inverse variance of area effects or 
Γ(1,  0.1) when only unstructured effects are used (see Figure <ref>). 


In contrast, design priors are obtained from the marginal posterior density of θ obtained from fitting ℳ on a simulated pilot sample s_0 with sampling fraction f=0.01. Where prior knowledge expressed in fitting priors is very informative, the sampling fraction can be relatively small and still produce relatively narrow design priors.





 §.§ Design-based simulation
 
A second series of simulations is carried out with a design-based procedure: samples are taken directly from the census tables using a multinomial law. Unlike with the model-based simulation, results no longer assume that the models are true. Doing so makes it possible to verify the validity of the sample size obtained under the previously described model-based SSD procedure. We also take the opportunity to estimate the relative efficiency in a series of scenarios of ascending complexity, from the simple design-based estimator to an HB predictor borrowing strength from a covariate and an explicit spatial covariance structure. These scenarios notated  S1 to S4 are summarised in Table <ref>. To this end, we simulate sampling and estimating each one and compare it with the true population value. 




One hundred simple random samples s_b are drawn from U  with simulated sizes  n_jdb∼Binomial(N_jdb, f) where f is set by the experiment, with no guarantee that all n_jdb >0. We simulate sample headcounts such that:

    y_jdb∼Binomial(n_jdb, N_jdb^-1 Y_jdb)


Samples s_b are used to produce model-based estimates Y_jdb of prevalence Y_jdb =   Y_jdb/N_jdb of LLTI every area making up the population as a function of characteristics of this population.

To reduce the computational burden of this simulation study, we only estimate relative efficiency under two sampling fractions: 2 and 4 per cent. The number of iterations is set to B=400 for each combination of a scenario and sampling fraction, totalling 2,400 procedures of model estimation and prediction. For scenarios S2–S5, and for each of the chosen sampling fractions f, we compute measures of accuracy:

	
  * the RMSE (RMSE) 
	
    RMSE_jd = √(B^-1∑_b=1^B(Y_jdb- Y_jdb)^2)

	
  * the bias
	
    Bias_jd = B^-1∑_b=1^B(Y_jdb- Y_jdb)

	
  * the absolute relative bias (ARB)
	
    ARB_jd = B^-1∑_b=1^B|Y_jdb- Y_jdb|Y_jdb

	
  * the relative RMSE or RSE (RSE) 
	
    RSE_jd = RMSE_jdY_jd

	
  * the RSE's relative bias (RSEB) (where postMSE is the posterior variance)
	
    RSEB_jd = B^-1∑_b=1^BRSE_jdb-  RSE_jdbRSE_jdb
       RSE_jd = √(postMSE_jd)Ŷ̅̂_jd



These measures are mapped and averaged across small areas d to be reported in tables. For scenario S1, the bias of sample means is by definition zero and we calculate the MSE and RSE using the variance formula of the sample proportion:

    Var(y̅_jd) = Y(1-Y)n_jd =  Y(1-Y) f  N_jd







§ RESULTS
 


 §.§ SSD procedure

The basic model (<ref>) is fitted to pilot sample s_0 of ESS n=146,574, with 1,952 cells jd containing no observations. This means no data is available for 1.3 per cent of the total number of cells, or 8.8 per cent of cells where the minimum quality standard applies. The fitted model achieves an acceptable DIC of 33,787. The SSD procedure summarised in Table <ref> produces an effective sample fraction requirement of 2.6875 per cent, equivalent to an ESS of 394,103. Assuming the procedure is correct, this solution is considered to fall at the most 0.1875 percentage points from the true effective sampling fraction achieving a negligible risk of cell suppression. Further iterations would reduce this distance. The algorithm was guided by an estimator of the loss function ℓ which depends on predictors θ to determine which cells should be checked for the quality standard. Any bias on θ can potentially increase or decrease the suppression risk to a substantial degree. In this case study, we note that the estimator ℓ̂ fails to detect a risk of suppression at step k=3, which leads the algorithm to select a lower interval of sample fractions that it would have done if it had been driven by actual loss ℓ  calculated using the census figures Y_jd.



Traditional survey sampling methods can determine the most cost-effective strategy to achieve this target ESS. A basic household survey like the EHS is based on a multistage design sampling one household per dwelling per address selected in the address sampling frame, and interviewing all eligible members (by proxy if necessary). Using the 2010/11 microdata sample for England of 11,188 effective responses, we estimated the DEFF to 1.16 for the LLTI question. This implies that survey planners should aim for a final sample size at least 16 per cent higher than the ESS found above, to account for the design effect. Yet many surveys have additional sampling stages which can increase the DEFF, for instance selecting postcode sectors as primary sampling units and then addresses/dwellings/households as secondary sampling units. This is for instance the case of the Health Survey for England for which we estimated a DEFF of 1.44 using a microdata sample from 2011 <cit.>. The DEFF can easily be used in area-level SAE models in a variance function, or in unit-level generalised linear mixed models with the ESS method described by <cit.>.



 §.§ Design-based simulation

We verify the validity of the small area model used in this case study as well as the reliability of the loss function estimator across iterations of the algorithm. Design-based simulation results are presented in Tables <ref> to <ref>. 

Results provide evidence that the relative posterior variances obtained from the spatial model are inflated by 20 to 40 per cent in basic scenario S2 (Table <ref>) and by as much as 40 to 50 per cent in scenario S3 (Table <ref>). This is likely to make the SSD procedure very conservative as many cells will be unnecessarily suppressed. This is a likely sign of poor model specification. Simulations for other scenarios are reported here to quantify the relative efficiency of the different model- and design-based methods. Efficiency savings from using a particularly model or from an increase in sample size can be derived as the ratio of RMSE of two scenarios. Comparing the scenarios presented in Table <ref>, we find that even a basic model without area covariates (see S1–S2, Tables <ref> and <ref>) achieves a reduction in RMSE by about 49 per cent overall, and up to 74 per cent for youngest age groups (f=0.02). This comes at a very reasonable computational cost and with a bias of between 0.1 and 0.7 percentage point. With a much larger sample (f=0.04), the efficiency gain comes down to approximately 35 per cent.

The addition of covariates delivers strong efficiency gains (see S3, Tables <ref>) with an reduction in RMSE by 72  (f=0.02) and 62 per cent (f=.04). This includes a reduction in bias from 0.02 to below 0.01 percentage point. 

Comparing results for f=0.02 and f=0.04, we note that doubling the sample size has diminishing returns as the methods becomes more complicated. Whilst it reduces the standard deviation of design-based estimates S1 by 29 per cent, the RMSE of estimates S2 is only reduced by 12 per cent. As for estimates S3 and S4 the reduction by less than 3 and 4 per cent respectively is not material for an augmentation of this magnitude. This shows that RMSE as a function of sample size is already relatively horizontal in the regions of sample sizes we are investigating, which is not the case for design-based estimators. Regardless of the sample size, bias remains high. Although the average bias
across areas remains of the order of 0.1 to 0.3 percentage points, the more meaningful ARB metric reveals an average absolute bias of 5 to 20 per cent of the target parameter. This represents almost all the total RSE. This bias is not very sensitive to the sampling fraction or the type of between-area variance structure used.

The working model was designed to borrow strength by assuming that areas situated near each other were more similar. We thus examined simulation results to verify whether this assumption holds everywhere equally. Figures <ref> and <ref> present the RMSE and bias for age-specific MSOA prevalence proportion computed from S3 with f=0.020. Hotpots of high RMSE and bias of 1 to 2 percentage points are visible which exhibit a spatial pattern. This could signal outlier areas, or areas in which that covariates are poorly measured and making the linear predictor particularly biased. The convolution prior used here may not fully and authentically reproduce the spatial pattern of the health status under consideration.

















§ FINAL REMARKS
 
This paper has examined the feasibility of a design approach to SAE model-based estimation, by proposing a fully Bayesian design algorithm taking a working model and a pilot sample as inputs. The decision rule is based on a new loss function of the sampling fraction: the weighted cell suppression rate. This approach is intended for survey planners and public health professionals to consider an SAE problem at the survey design phase, rather than using it strictly for post hoc analyses. The availability of historical survey data is likely to allow survey planners to implement similar approaches at a small cost. The Bayesian decision-analytic approach is well suited to make use of whatever information is available—in this case, a pilot sample—but it is perfectly possible to form a comparable design prior on the basis of good information on model random component (random effect variances, spatial fraction). The SSD algorithm then operates on a simulated population that is only as realistic and informative as the design prior is.

Our approach is a direct continuation of previous model-based design decisions considered by <cit.>, <cit.> and <cit.>. It establishes a new link to official statistics requirements with a new loss function directly relevant to quality criteria present in most statistical codes of practice internationally. This loss function can be calculated for different social subgroups in situations where estimates of equal precision are needed for minority populations, particularly when equality impact assessments demand equal statistical precision for official statistics across different communities. 

The implementation illustrated in our case study greatly simplifies the design work by expressing the output in terms of effective sampling fraction or ESS, which integrates well with existing survey design tools, knowledge around the DEFF of various sampling designs, and work to address sample attrition caused by nonresponse.  The case study is designed around random sampling stratified by small area with allocation proportional to size, but more sophisticated approaches can be considered in relation to optimal unequal probability sampling using inferential priority weights as proposed by <cit.>, incorporating on known population characteristics across areas (for instance ethnic diversity or age characteristics). It is also relevant to other model-based estimation techniques such as generalised structure preserving estimation <cit.>.

While the computational burden of this method can remain high in studies involving a large number of cells, the combination of INLA and the binary search logic brings these costs to levels that are no longer absurd with modern equipment and high performance cloud computing. Most importantly, the binary search allows planners to anticipate the exact number of iterations needed to reach a prespecified solution precision. They can therefore budget the algorithm's computational requirements.

There are limitations to this approach which make it highly dependent on model misspecification and modelling priors which, together with model selection, arguably constitute the three main challenges to model-based estimation. Given known difficulties in identifying spatial models, it is clear that such a procedure is not likely to be reliable outside of models well established across years of data collection and model testing. While this is not unrealistic with government surveys, this does represent a non-negligible design cost which is difficult to incorporate explicitly in planning.

Our case study conducted on a real population also illustrates that the computation of the loss function being dependent on the reliability of the model, model diagnostics and validation are essential to the reliability and stability of the procedure. Results from the design-based simulation studies show the influence of typical model misspecification on estimation bias and how it clusters in space or social groups. Prior formation also remains, as in many other areas of statistics, a challenge. Recent work in this area by <cit.> and <cit.> expresses some of the important obstacles faced by practitioners in this area.

This work tends to support other research showing the feasibility of reconciling Bayesian inference and survey sampling, the intersection of which, in this case, is expressed in terms of ESS. This also opens perspectives to combine survey sampling with some very strong developments witnessed in recent years in the design of experiments literature, particularly around adaptive trial designs.
 



§ ACKNOWLEDGEMENTS

The author would like to thank Professor Sujit Sahu and Doctor Marta Blangiardo for their helpful comments on an early draft. The author acknowledges the use of the IRIDIS High Performance Computing Facility, and associated support services at the University of Southampton, in the completion of this work.


§ FUNDING

This work was supported by an Economic & Social Research Council, Advanced Quantitative Methods doctoral studentship [reference number 1223155].


§ ABBREVIATIONS


 
	
	
	AICAkaike information criterion
	AR1first-order autoregressive
	ARBabsolute relative bias
	AUCarea under the curve
	BMIbody mass index
	CARconditional autoregressive
	CCGClinical Commissioning Group
	CEcommunal establishment
	CSEWCrime Survey for England and Wales
	CVcoefficient of variation
	DEFFdesign effect
	DHSDemographic and Health Survey
	DICdeviance information criterion
	EBempirical best
	EBLUPempirical best linear unbiased predictor
	EHSEnglish Housing Survey
	EPPempirical plug-in predictor
	EPSEMequal probability of sampling method
	ESSeffective sample size
	FRSFamily Resources Survey
	GISgeographical information system
	HBhierarchical Bayes
	HSCICHealth & Social Care Information Centre 
	HWBHealth and Wellbeing Board
	ICARintrinsic conditional autoregressive
	ICDInternational Classification of Diseases 
	IHSIntegrated Household Survey
	INLAintegrated nested Laplace approximation
	ISARindirectly standardised emergency hospital admissions rate
	JSNAjoint strategic needs assessment
	LADlocal authority district
	LFSLabour Force Survey
	LISAlocal indicator of autocorrelation
	LLBLeroux-Lei-Breslow 
	LLTIlimiting long-term illness
	LOSLife Opportunities Survey
	MORmedian odds ratio
	MRPmultilevel regression and poststratification
	MSEmean squared error
	MSOAmiddle layer super output area
	MYPEmid-year population estimate
	NEXnormal exchangeable 
	NHSNational Health Service
	ONSUK Office for National Statistics
	PCTPrimary Care Trust
	QLFSQuarterly Labour Force Survey
	Q–Qquantile-quantile
	RBrelative bias
	RMSEroot mean squared error
	ROCreceiver operating characteristics
	RSErelative standard error
	RW1first-order random walk
	SAEsmall area estimation
	SARsimultaneous autoregressive
	SILCSurvey on Income and Living Conditions
	SMRdirectly age-standardised mortality rate
	SRHself-rated health
	SSDsample size determination
	SYNsynthetic predictor
	UPSEMunequal probability of sampling method






