Estimation of domain discontinuities using Hierarchical Bayesian Fay-Herriot models
Section 3. Methods

• Previous
• Next

3.1   Small area estimation for domain discontinuities

Testing hypotheses about differences between estimates of a finite population parameter observed under different survey processes implies the existence of measurement errors. Therefore a measurement error model is required to explain systematic differences between survey estimates for the same population parameter observed under two different survey approaches. Let ${\theta }_i$ denote the population parameter of domain $i\mathrm{<mo>=</mo>}\mathrm{1,}\dots \mathrm{,}m.$ Let $y_i^r$ and $y_i^a$ denote the observed value for ${\theta }_i$ in the case of a complete enumeration under the regular approach and alternative approach, respectively. Direct estimates for $y_i^r$ and $y_i^a$ are obtained with the GREG estimator based on the samples assigned to the regular and alternative survey and are denoted as ${\widehat{y}}_i^r$ and ${\widehat{y}}_i^a$ respectively.

The relation between the observed values under a complete enumeration and the real population parameter is:   $y_i^q\mathrm{<mo>=</mo>}{\theta }_i+{\lambda }_i^q,$ $i\mathrm{<mo>=</mo>}\mathrm{1,}\dots \mathrm{,}m,$ $q\mathrm{<mo>=</mo>}r\mathrm{,}a,$ with ${\lambda }_i^q$ the real measurement bias if ${\theta }_i$ is measured with survey approach $q.$ Without any external information, it is not possible to estimate ${\lambda }_i^q.$ From the sample data only the relative bias, say ${\Delta }_i\mathrm{<mo>=</mo>}{y}_i^r-y_i^a\mathrm{<mo>=</mo>}{\lambda }_i^r-{\lambda }_i^a$ is identifiable. Direct estimates for these discontinuities are obtained from the survey data as the contrast between the GREG estimates, i.e., ${\widehat{\Delta }}_i\mathrm{<mo>=</mo>}{\widehat{y}}_i^r-{\widehat{y}}_i^a.$

In the case of the Dutch CVS the sample size of the regular survey is large enough to obtain sufficiently precise direct estimates for the planned domains, since the sample is designed to publish official statistics for these domains. The sample assigned to the alternative survey for the parallel run has only a size of one third of the regular sample size, which is insufficient to obtain precise direct estimates for the planned domains. In an earlier paper (van den Brakel et al., 2016) univariate FH models were developed to obtain more precise predictions for the domain parameters observed with the small sample size assigned to the alternative survey approach using auxiliary variables derived from three different sources. The first source contains demographic variables derived from the Municipal Basis Administration (MBA), which is an administration of all people residing in the Netherlands. The second source contains related variables available in the Police Register of Reported Offences (PRRO). The third source, which is unique in the case of a parallel run, contains direct estimates for the same variables observed under the regular survey, which are sufficiently precise at least for the planned domains like police districts. The direct estimates from the regular survey are often selected as auxiliary variables for these univariate FH models. This comes not as a surprise since these are survey estimates for the same population parameters. Although measured with a different survey process, strong positive correlations can be expected. Strong improvements of the precision of small domain prediction are indeed found if the set of potential auxiliary variables, i.e., from MBA and PRRO, is extended with the direct estimates from the regular CVS.

In this application the sampling error in the auxiliary variables that come from the regular CVS can be ignored in the FH model, since the sample size and therefore the sampling error for these domains is more or less equal for the domains (Ybarra and Lohr, 2008). This implies that the variance component of the random domain effects is inflated with the sampling error of the auxiliary variables, which is fine as long as the sampling error does not differ between domains. In most applications this is not the case and the methods proposed by Ybarra and Lohr (2008) should be used to account for sampling error in the auxiliary variables.

FH multilevel models can be fitted under a frequentist approach using EBLUP or under an HB approach (Rao and Molina, 2015). In van den Brakel et al. (2016) the HB approach is preferred over the EBLUP, since the strong auxiliary information in the fixed effect part of the model often results in zero estimates for the variance component of the random domain effects, giving too much weight to the synthetic regression part and too little weight to the direct estimates in the EBLUP, (Bell, 1999; Rao and Molina, 2015). This problem can also be overcome with adjusted maximum likelihood estimation, see e.g., Li and Lahiri (2010) and Hirose and Lahiri (2018).

Let ${\tilde{y}}_i^a$ denote the HB prediction for domain $i$ under the alternative approach. Now domain discontinuites are obtained by ${\tilde{\Delta }}_i\mathrm{<mo>=</mo>}{\widehat{y}}_i^r-{\tilde{y}}_i^a.$ Using direct estimates of the regular survey as auxiliary variables in the fixed part of the FH model for the alternative survey considerably increases the complexity of the variance estimation for the discontinuities. The variance of ${\tilde{\Delta }}_i$ can be expressed as $\text{Var}\left({\tilde{\Delta }}_i\right)\mathrm{<mo>=</mo>}\text{Var}\left({\widehat{y}}_i^r\right)+\text{MSE}\left({\tilde{y}}_i^a\right)-2\text{cov}\left({\widehat{y}}_i^r\mathrm{,}{\tilde{y}}_i^a\right).$ The use of ${\widehat{y}}_i^r$ or related sample estimates as auxiliary variables to predict ${\tilde{y}}_i^a,$ results in non-zero values for $\text{cov}\left({\widehat{y}}_i^r\mathrm{,}{\tilde{y}}_i^a\right)$ that cannot be ignored. To approximate $\text{Var}\left({\tilde{\Delta }}_i\right),$ van den Brakel et al. (2016) proposed an approximately design-unbiased estimator for $\text{cov}\left({\widehat{y}}_i^r\mathrm{,}{\tilde{y}}_i^a\right)$ and $\text{Var}\left({\widehat{y}}_i^r\right),$ while the $\text{MSE}\left({\tilde{y}}_i^a\right)$ is approximated with the posterior variance of the HB domain predictions. A major disadvantage of this approach is that model-based and design-based uncertainty measures are intertwined. On the one hand, the MSE’s for ${\tilde{y}}_i^a$ are approximated with their posterior variances. On the other hand, the covariances between ${\widehat{y}}_i^r$ and ${\tilde{y}}_i^a$ are approximated from a design-based perspective. Consequently, naive application of this approach may give negative variance estimates for the discontinuities. This drawback has been solved using a design-based approximation for the $\text{MSE}\left({\tilde{y}}_i^a\right),$ resulting in a fully design-based approximation for the uncertainty of the estimated domain discontinuities.

In this paper a full HB framework for estimating domain discontinuities is proposed as an alternative by developing a bivariate FH model for the domain parameters observed under both the regular and alternative approach. The advantage of this approach is that it improves the precision of both the direct estimates of the regular and alternative domain estimates by borrowing strength from other domains and both surveys. Negative variance estimates for the estimated domain discontinuities are precluded by definition under this multivariate HB framework. This method is compared with a simple alternative, namely a univariate FH model for the direct estimates of the discontinuities. As mentioned in the introduction, it is anticipated that it might be hard to find covariates in the available registers that explain the discontinuities. An advantage of both models is that they avoid the complications of accounting for sampling error in the auxiliary variables, which is necessary if the survey estimates of the regular survey are used as covariates in univariate FH models and the sampling error differs between domains.

3.2   Bivariate Fay-Herriot model

A bivariate version of the FH model (Fay and Herriot, 1979) starts with a measurement model for the two GREG estimates observed in each domain:

$${\widehat{y}}_i\mathrm{<mo>=</mo>}{y}_i+e_i\mathrm{,}i\mathrm{<mo>=</mo>}\mathrm{1,}\dots \mathrm{,}m\mathrm{,}(3.1)$$

with $y_i\mathrm{<mo>=</mo>}{\left(y_i^r\mathrm{,}{y}_i^a\right)}^t,$ ${\widehat{y}}_i$ a vector containing the GREG estimates of $y_i$ and $e_i\mathrm{<mo>=</mo>}{\left(e_i^r\mathrm{,}{e}_i^a\right)}^t$ a vector with the sampling errors of ${\widehat{y}}_i$ for which it is assumed that

$$e_i\stackrel{\text{ind}}{{\displaystyle ~}}N\left(0_2\mathrm{,}{\Psi }_i\right),i\mathrm{<mo>=</mo>}\mathrm{1,}\dots \mathrm{,}m\mathrm{.}(3.2)$$

Here $0_2$ is a 2 dimensional column vector with each element equal to zero. Since the sample for the regular and alternative survey are drawn independently, it is assumed that ${\Psi }_i\mathrm{<mo>=</mo>}\text{Diag}\left({\psi }_i^r\mathrm{,}{\psi }_i^a\right)$ where ${\psi }_i^q$ is the design variance of ${\widehat{y}}_i^q.$ It is also assumed that these design variances are known although they are replaced by their estimates in practice. The true domain parameters are modelled with a multilevel model. For the fixed effects it is assumed that the regular and alternative approach share the same covariates. In the most general case the regression coefficients for the fixed part are different for both variables i.e., $y_i^q\mathrm{<mo>=</mo>}{x}_i^t{\beta }^q+{\nu }_i^q,$ with $x_i\mathrm{<mo>=</mo>}{\left(x_{i1}\mathrm{,}\dots \mathrm{,}{x}_{ip}\right)}^t$ a $p$ -vector with covariates of domain $i,$ ${\beta }^q$ a $p$ -vector of regression coefficients, which are equal over the domains but might be different between the two survey approaches. It is assumed that $x_{i1}\mathrm{<mo>=</mo>}1$ corresponds to the intercept. Furthermore ${\nu }_i^q$ are random domain effects. This gives rise to the following bivariate multilevel model for the two domain parameters:

$$y_i\mathrm{<mo>=</mo>}{X}_i\beta +{\nu }_i\mathrm{,}i\mathrm{<mo>=</mo>}\mathrm{1,}\dots \mathrm{,}m\mathrm{,}(3.3)$$

where $X_i\mathrm{<mo>=</mo>}{I}_2\otimes x_i^t,$ $\beta \mathrm{<mo>=</mo>}{\left({\beta }^{r^t},{\beta }^{a^t}\right)}^t,$ $I_2$ a 2 dimensional identity matrix, and ${\nu }_i\mathrm{<mo>=</mo>}{\left({\nu }_i^r\mathrm{,}{\nu }_i^a\right)}^t.$ For the random domain effects it is assumed that

$${\nu }_i\stackrel{\text{IID}}{{\displaystyle ~}}N\left(0_2\mathrm{,}\Sigma \right)\mathrm{,}i\mathrm{<mo>=</mo>}\mathrm{1,}\dots \mathrm{,}m\mathrm{,}(3.4)$$

with $\Sigma$ a general $2\times 2$ covariance matrix for the random domain effects. Inserting (3.3) into (3.1) gives:

$${\widehat{y}}_i\mathrm{<mo>=</mo>}{X}_i\beta +{\nu }_i+e_i\mathrm{,}i\mathrm{<mo>=</mo>}\mathrm{1,}\dots \mathrm{,}m\mathrm{,}(3.5)$$

with model assumptions (3.2) and (3.4).

Since the number of domains in this application is small, it is important to select parsimonious models. One way to reduce model complexity is to assume that the regression coefficients are equal for both survey approaches. In this case a dummy indicator, say ${\delta }_i,$ is introduced, which is equal to zero for the regular survey and equal to one for the alternative survey. In this case $x_{i1}\mathrm{<mo>=</mo>}1$ corresponds to the overall intercept and $x_{i2}\mathrm{<mo>=</mo>}{\delta }_i$ is the indicator whose coefficient measures the differences between intercepts of the variables observed under both surveys. So $y_i^q\mathrm{<mo>=</mo>}{x}_i^t\beta +{\nu }_i^q$ and in (3.3) $X_i\mathrm{<mo>=</mo>}{x}_i^t,$ and $\beta$ a vector with the corresponding regression coefficients. As a result, two versions for the fixed effects are considered: 

• FE_uq: A fixed effect model where the regular and alternative approach share the same covariates, but have different regression coefficients. In this case, domain discontinuities are given by

$${\Delta }_i\mathrm{<mo>=</mo>}{\displaystyle \sum _{j\mathrm{<mo>=</mo>\; 1}}^p{x}_{i\mathrm{,}j}\left({\beta }_j^r-{\beta }_j^a\right)}+\left({\nu }_i^r-{\nu }_i^a\right)\mathrm{.}(3.6)$$

• FE_eq: A more parsimonious version for the fixed effect component by assuming that the regression coefficients are equal for the regular and alternative approach. In this case domain discontinuities are given by

$${\Delta }_i\mathrm{<mo>=</mo>}-{\beta }_2+\left({\nu }_i^r-{\nu }_i^a\right)\mathrm{,}(3.7)$$

• with ${\beta }_2$ the regression coefficient for $x_{i2}\mathrm{<mo>=</mo>}{\delta }_i.$

The following covariance structures for the random domain effects are considered: 

• RE_f: A full covariance matrix $\Sigma$ for the random domain effects. Positive correlation between the random domain effects will further increase the precision of the estimates for the domain discontinuities since domain estimates borrow strength not only from different domains but also across the two surveys.
• RE_d: A diagonal covariance matrix with separate variances for the regular and alternative approach, i.e.: $\Sigma \mathrm{<mo>=</mo>}\text{Diag}\left({\sigma }_r^2\mathrm{,}{\sigma }_a^2\right).$ This covariance structure in combination with model FE_uq comes down to applying a univariate FH model to both surveys separately. In this case models only use sample information from other domains within the same survey but not across the two surveys to improve the precision of the estimates for domain discontinuities.
• RE_s: A diagonal covariance matrix with equal variances for the regular and alternative approach, i.e.: $\Sigma \mathrm{<mo>=</mo>}{\sigma }^2{I}_2.$

3.3   Univariate Fay-Herriot model for domain discontinuities

The univariate FH model for the direct estimates of the discontinuities starts with defining a measurement error model for the GREG estimates of the discontinuities:

$${\widehat{\Delta }}_i\mathrm{<mo>=</mo>}{\Delta }_i+z_i(3.8)$$

with ${\Delta }_i\mathrm{<mo>=</mo>}{y}_i^r-y_i^a$ the true discontinuity of domain $i$ under a complete enumeration of the population under both approaches, ${\widehat{\Delta }}_i\mathrm{<mo>=</mo>}{\widehat{y}}_i^r-{\widehat{y}}_i^a$ the GREG estimate for ${\Delta }_i$ based on the parallel run and $z_i\mathrm{<mo>=</mo>}{e}_i^r-e_i^a$ the sampling error of ${\widehat{\Delta }}_i.$ It is assumed that $z_i$ $\simeq$ $N\left(\mathrm{0,}{\psi }_i^r+{\psi }_i^a\right)$ and that the design variances of the sampling errors are known. For ${\Delta }_i$ the following linear model is assumed:

$${\Delta }_i\mathrm{<mo>=</mo>}{x}_i^t\beta +v_i\mathrm{,}(3.9)$$

with $v_i$ $\simeq$ $N\left(\mathrm{0,}{\sigma }_v^2\right).$ Inserting (3.9) into (3.10) gives:

$${\widehat{\Delta }}_i\mathrm{<mo>=</mo>}{x}_i^t\beta +v_i+z_i.(3.10)$$

3.4   Estimation of the bivariate Fay-Herriot model

The models developed in Subsections 3.2 and 3.3 are fitted with a HB approach using Markov Chain Monte Carlo (MCMC) sampling. In particular the Gibbs sampler is used. The following priors are used for the model parameters and hyperparameters. For the regression coefficients uniform improper priors are assumed, i.e., $\beta ~1.$ For the random domain effects a multivariate normal prior is used: $\nu |\Sigma$ $\simeq$ $N\left(0_{2m}\mathrm{,}{I}_m\otimes \Sigma \right).$

In the case of a full covariance matrix for the random domain effects, the prior for $\Sigma$ is taken to be a scaled inverse Wishart distribution (O’Malley and Zaslavsky, 2008). This distribution is obtained by writing $\Sigma \mathrm{<mo>=</mo>}\text{Diag}\left(\xi \right)\tilde{\Sigma }\text{Diag}\left(\xi \right),$ with $\xi \mathrm{<mo>=</mo>}{\left({\xi }^r\mathrm{,}{\xi }^a\right)}^t$ and assuming a standard normal distribution for ${\xi }^r$ and ${\xi }^a,$ i.e., ${\xi }^x$ $\simeq$ $N\left(\mathrm{0,}1\right),$ $\left(x\in \left(a\mathrm{,}r\right)\right)$ and an inverse Wishart distribution for $\tilde{\Sigma },$ i.e., $\tilde{\Sigma }$ $\simeq$ $\text{Inv}-\text{Wish}\left(v_v\mathrm{,}{\Phi }_v\right),$ with $v_v\mathrm{<mo>=</mo>}d+1$ degrees of freedom, with $d$ the dimension of $\tilde{\Sigma }$ which is equal to 2 in this application, and scale parameter ${\Phi }_v\mathrm{<mo>=</mo>}{I}_2.$ In the case of a diagonal covariance matrix for the random domain effects, the priors for ${\sigma }_q$ (in the case of unequal variances) or $\sigma$ (in the case of equal variances) are half-Cauchy distributions. These are more robust prior distributions than the more commonly used inverse chi-squared distribution (Gelman, 2006). The inverse chi-squared distribution might be informative, even in the case of small scale and shape parameters. In addition convergence problems might occur with the Gibbs sampler. Both problems are largely avoided with a redundant multiplicative parametrization of the random effects (Gelman, 2006; Gelman, Van Dyk, Huang and Boscardin, 2008; Polson and Scott, 2012). See van den Brakel and Boonstra (2018) for more details on the priors of this model.

Let $\widehat{y}$ denote the $2m$ column vector obtained by stacking the $m$ column vectors ${\widehat{y}}_i,$ and $X$ the matrix obtained by stacking the matrices $X_i.$ In the case of unequal regression coefficients, $X$ is a $2m\times 2p$ matrix. In the case of equal regression coefficients, $X$ is a $2m\times p$ matrix. The likelihood function can be written as

$$p\left(\widehat{y}|\theta \right)\mathrm{<mo>=</mo>}N\left(X\beta +\nu \mathrm{,}\Psi \right)\mathrm{,}(3.11)$$

with $\Psi \mathrm{<mo>=</mo>}{\oplus }_{i\mathrm{<mo>=</mo>\; 1}}^m{\Psi }_i$ a $2m\times 2m$ diagonal matrix with the design variances of the direct estimates $\widehat{y}$ and $\theta$ a vector containing all model parameters. The joint prior distribution $p\left(\text{}\theta \text{}\right)$ equals the product of the aforementioned priors. The posterior distribution of $\theta$ is proportional to the joint density, i.e., $p\left(\theta |\widehat{y}\right)\propto p\left(\text{}\theta \text{}\right)p\left(\widehat{y}|\theta \right).$ The model is fitted using the Gibbs sampler Geman and Geman (1984); Gelfand and Smith (1990). The full conditional distributions used in the Gibbs sampler are specified in van den Brakel and Boonstra (2018).

For each model considered, the Gibbs sampler is run in three independent chains with randomly generated starting values. The length of each chain after the burn-in period for each run is 10,000 iterations. This gives 30,000 draws to compute estimates and standard errors. The convergence of the MCMC simulation is assessed using trace and autocorrelation plots as well as the Gelman-Rubin potential scale reduction factor (Gelman and Rubin, 1992), which diagnoses the mixing of the chains. The diagnostics suggest that all chains converge well within 500 draws. The estimated Monte Carlo simulation errors are small compared to the posterior standard errors for all parameters, so that the number of draws are more than sufficient for our purposes.

The estimands of interest are expressed as functions of the parameters, and applying these functions to the MCMC output for the parameters results in draws from the posteriors for these estimands. Domain predictions for the target variables under the bivariate FH model are obtained as the posterior means approximated by the Gibbs sampler output and are denoted as ${\tilde{y}}_i^{q\mathrm{,}b\text{FH}}.$ Domain predictions for the discontinuities are obtained as the posterior means of (3.6) or (3.7) approximated by the Gibbs sampler output and are denoted as ${\tilde{\Delta }}_i^{b\text{FH}}.$ Mean squared errors for ${\tilde{y}}_i^{q\mathrm{,}b\text{FH}}$ and ${\tilde{\Delta }}_i^{b\text{FH}}$ are obtained as posterior variances approximated from the Gibbs sampler output.

The methods are implemented in R using the mcmcsae R-package (Boonstra, 2020).

3.5   Pooling design variances

Estimates for the design variances ${\psi }_i^r$ and ${\psi }_i^a$ are available from the GREG estimator and are used as if the true design variances are known. This is a standard assumption in small area estimation. Therefore it is important to provide reliable estimates for these design variances. For the regular survey the variance estimates of the GREG estimates are considered to be reliable enough to be used in the FH model. For the alternative survey the estimates of the design variances are unreliable and therefore smoothed to improve their stability of the estimates of ${\psi }_i^a.$ Under the assumption that the population variances of the GREG residuals under the alternative approach are equal accross domains, the analysis-of-variance type of pooled variance estimator is used:

$${\psi }_i^a\mathrm{<mo>=</mo>}\frac{1-f_i^a}{n_i^a}\frac{1}{n^a-m}{\displaystyle \sum _{i=1}^m\left(n_i^a-1\right)S_{i;\text{GREG}}^{a^2}}$$

with $f_i^a$ the sample fraction in domain $i$ of the alternative survey, $n_i^a$ the number of respondents in domain $i$ under the alternative survey, $n^a\mathrm{<mo>=</mo>}{\displaystyle {\sum }_{i\mathrm{<mo>=</mo>\; 1}}^m{n}_i^a},$ and $S_{i\mathrm{<mo>;</mo>}\text{GREG}}^{a^2}$ the estimated population variance of the GREG residuals.

Alternatively, variance estimates of the direct estimates can be smoothed by modeling the variance estimates along with the GREG estimates themselves, (You and Chapman, 2006) and Sugasawa, Tamae and Kubokawa (2017). Their approach can be traced back to Arora and Lahiri (1997). Another possibility is to smooth variance estimates by applying generalized variance functions, (Wolter (2007), Chapter 7, and Hawala and Lahiri (2018)).

3.6   Model selection and evaluation

Frequently applied model selection criteria in HB settings are the Widely Applicable Information Criterion or Watanabe-Akaike Information Criteria (WAIC) (Watanabe, 2010, 2013) and the Deviance Information Criteria (DIC) (Spiegelhalter, Best, Carlin and van der Linde, 2002). They are popular because they are easy to compute from MCMC simulation output and because of their ability to make a reasonable tradeoff between model fit and model complexity. The WAIC is seen as an improvement on the DIC since the latter can produce negative estimates for the effective number of parameters and it is not defined for singular models (Vehtari, Gelman and Gabry, 2017). The penalty used for model complexity in DIC and WAIC is closely related to the effective number of parameters proposed by Hodges and Sargent (2001) for linear multilevel models where each fixed effect contributes one degree of freedom and the random effects contribute a value in the range between zero and $m,$ depending on the size of the variance component. As follows from the definition of WAIC, models with lower WAIC values are preferred. The WAIC estimates are uncertain and an approximation of its standard error is provided by Vehtari et al. (2017) equation (23) and can be computed using R package $loo$ (Vehtari, Gelman and Gabry, 2015).

Covariates are selected from the set of auxiliary variables listed in van den Brakel and Boonstra (2018) using a step-forward selection procedure. Various models are compared using the aforementioned WAIC estimates. From the set of potential covariates, the covariate with the lowest WAIC value is selected in the model. This selection process is iteratively repeated as long as adding a new covariate further decreases the WAIC value. In this application, this step-forward selection procedure, further abbreviated as step-WAIC, often results in models with a large number of covariates. Since the WAIC values are estimates that contain error, it appears that it might not be desirable to minimize the WAIC by adding covariates to the model as long as it reduces the point estimates of the WAIC. As an alternative we applied a step-forward selection procedure where covariates are added to the model as long as a new covariate decreases the WAIC with a value that exceeds the estimated standard error of the WAIC. This method will be referred to as step-WAIC-se.

The step-forward selection procedure is applied to each of the six different combinations of the two fixed effect versions (FE_uq and FE_eq) and the three covariance structures of the random component (RE_f, RE_d, and RE_s). From the resulting six models the one with the lowest WAIC value is selected. Model adequacy of these six selected models is evaluated with posterior predictive checks. This implies that replicate data sets, simulated from the posterior predictive distribution are compared with the originally observed data to study systematic discrepancies and to evaluate how well the selected model fits the observed data (Gelman, Carlin, Stern, Dunson, Vehtari and Rubin, 2004). Posterior predictive $p$-values are calculated for six different tests that evaluate particular aspects of the posterior predictive distribution. Posterior predictive $p$ -values for the domain discontinuities are defined as $p\mathrm{<mo>=</mo>}P\left(T\left({\widehat{\Delta }}^{\text{sim}}\mathrm{,}\Delta \right)\ge T\left(\widehat{\Delta }\mathrm{,}\Delta \right)|\widehat{y}\right),$ where ${\widehat{\Delta }}^{\text{sim}}\mathrm{<mo>=</mo>}{\left({\widehat{\Delta }}_1^{\text{sim}}\mathrm{,}\dots \mathrm{,}{\widehat{\Delta }}_m^{\text{sim}}\right)}^t$ are replicates of the observed discontinuities for the $m$ domains under the posterior predictive distribution, $\widehat{\Delta }\mathrm{<mo>=</mo>}{\left({\widehat{\Delta }}_1\mathrm{,}\dots \mathrm{,}{\widehat{\Delta }}_m\right)}^t$ the observed direct estimates for the $m$ domain discontinuities and $T({\widehat{\Delta }}^{\text{sim}},\Delta )$ (or $T(\widehat{\Delta }\mathrm{,}\Delta ))$ a test statistic that depends on ${\widehat{\Delta }}^{\text{sim}}$ (or $\widehat{\Delta })$ and unknown true values for the $m$ domain discontinuities $\Delta \mathrm{<mo>=</mo>}{\left({\Delta }_1\mathrm{,}\dots \mathrm{,}{\Delta }_m\right)}^t.$ Posterior predictive $p$ -values are estimated from the Gibbs sampler output as the average over the $S$ Monte Carlo samples

$$\widehat{p}\mathrm{<mo>=</mo>}\frac{1}{S}{\displaystyle \sum _{s\mathrm{<mo>=</mo>\; 1}}^{S}I\left(T\left({\widehat{\Delta }}^s\mathrm{,}{\Delta }^s\right)\ge T\left(\widehat{\Delta }\mathrm{,}{\Delta }^s\right)\right)\mathrm{,}}(3.12)$$

with $I\left(\text{}A\text{}\right)$ the indicator function with value one if the condition $A$ is fulfilled and zero otherwise, ${\widehat{\Delta }}^s\mathrm{<mo>=</mo>}{\left({\widehat{\Delta }}_1^s\mathrm{,}\dots \mathrm{,}{\widehat{\Delta }}_m^s\right)}^t$ the $m$ observed domain discontinuities in the $s^{\text{th}}$ replicate of the MCMC simulation and ${\Delta }^s\mathrm{<mo>=</mo>}{\left({\Delta }_1^s\mathrm{,}\dots \mathrm{,}{\Delta }_m^s\right)}^t$ the true values of the $m$ domain discontinuities in the $s^{\text{th}}$ replicate of the MCMC simulation. If a model fits the observed data adequately, then it is expected that $T\left(\widehat{\Delta }\mathrm{,}{\Delta }^s\right)$ is in the bulk of the histogram of the replicates $T\left({\widehat{\Delta }}^s\mathrm{,}{\Delta }^s\right).$ Therefore $p$ -values close to zero or one are indications of a poor fit with respect to that test statistic. In the expressions below, it is understood that $T_x,$ for $x\mathrm{<mo>=</mo>}\mathrm{1,}\dots \mathrm{,}6,$ is a function of $\left({\widehat{\Delta }}^s\mathrm{,}{\Delta }^s\right)$ or $\left(\widehat{\Delta }\mathrm{,}{\Delta }^s\right),$ depending on the component that is evaluated in (3.12). The following posterior predictive tests are defined (You, 2008):

1. A general goodness-of-fit test statistic $T_1\mathrm{<mo>=</mo>}{\displaystyle {\sum }_{i\mathrm{<mo>=</mo>\; 1}}^m{\left({\widehat{\Delta }}_i-{\Delta }_i\right)}^2/\text{Var}\left({\widehat{\Delta }}_i|{\Delta }_i\right)}.$ Here $\text{Var}\left({\widehat{\Delta }}_i|{\Delta }_i\right)\mathrm{<mo>=</mo>}{\psi }_i^r+{\psi }_i^a.$
2. $T_2\mathrm{<mo>=</mo>}\text{max}\left({\widehat{\Delta }}_i\right)$ and $T_3\mathrm{<mo>=</mo>}\text{min}\left({\widehat{\Delta }}_i\right)$ which are sensitive for deviations in the tails of the distribution.
3. $T_4\mathrm{<mo>=</mo>}\frac{1}{m}{\displaystyle {\sum }_{i\mathrm{<mo>=</mo>\; 1}}^m{\widehat{\Delta }}_i}\equiv \overline{\widehat{\Delta }},$ i.e., the mean which is sensitive for bias in the domain predictions.
4. $T_5\mathrm{<mo>=</mo>}\frac{1}{m-1}{{\displaystyle {\sum }_{i\mathrm{<mo>=</mo>\; 1}}^m\left({\widehat{\Delta }}_i-\overline{\widehat{\Delta }}\right)}}^2,$ i.e., the variance of the domain estimates, which is sensitive for e.g., overshrinkage.
5. $T_6\mathrm{<mo>=</mo>}\left|\text{max}\left({\widehat{\Delta }}_i\right)-\overline{\Delta }\right|-\left|\text{min}\left({\widehat{\Delta }}_i\right)-\overline{\Delta }\right|,$ with $\overline{\Delta }\mathrm{<mo>=</mo>}\frac{1}{m}{\displaystyle {\sum }_{i\mathrm{<mo>=</mo>\; 1}}^m{\Delta }_i}$ which is sensitive to asymmetry in the distribution.

• Previous
• Article main page
• Next