Modelling multiple hospital outcomes: the impact of small area and primary care practice variation

Background: Appropriate management of care – for example, avoiding unnecessary attendances at, or admissions to, hospital emergency units when they could be handled in primary care – is an important part of health strategy. However, some variations in these outcomes could be due to genuine variations in health need. This paper proposes a new method of explaining variations in hospital utilisation across small areas and the general practices (GPs) responsible for patient primary care. By controlling for the influence of true need on such variations, one may identify remaining sources of excess emergency attendances and admissions, both at area and practice level, that may be related to the quality, resourcing or organisation of care. The present paper accordingly develops a methodology that recognises the interplay between population mix factors (health need) and primary care factors (e.g. referral thresholds), that allows for unobserved influences on hospitalisation usage, and that also reflects interdependence between hospital outcomes. A case study considers relativities in attendance and admission rates at a North London hospital involving 149 small areas and 53 GP practices. Results: A fixed effects model shows variations in attendances and admissions are significantly related (positively) to area and practice need, and nursing home patients, and related (negatively) to primary care access and distance of patient homes from the hospital. Modelling the impact of known factors alone is not sufficient to produce a satisfactory fit to the observations, and random effects at area and practice level are needed to improve fit and account for overdispersion. Conclusion: The case study finds variation in attendance and admission rates across areas and practices after controlling for need, and remaining differences between practices may be attributable to referral behaviour unrelated to need, or to staffing, resourcing, and access issues. In managerial terms, the analysis points to the utility of formal statistical analysis of hospitalisation rates as a prelude to non-statistical investigation of primary care resourcing and organisation. For example, there may be implications for the location of staff involved in community management of chronic conditions; health managers may also investigate whether some practices have unusual populations (homeless, asylum seekers, students) that explain different hospital use patterns.


Background
Of particular importance in strategic management of the primary-acute care interface is identification of the sources of variation in hospital referrals and attendances, and detecting whether particular small areas and GP practices have above average emergency attendance and admission rates, especially where such variations are not related to acknowledged sources of health need. Containment of unplanned emergency admissions to, and attendances at, hospital emergency units is a major element in strategic management of the health demand generated by long term chronic illness, as hospital based (acute) care is relatively costly. If excess referral rates are indicated for particular practices this may have implications for clinical and organisational effectiveness in primary care -see, for example [1]. However, a statistical analysis of the kind undertaken here can only identify referral variations per se, and not necessarily make inferences on effectiveness.
In the recent past, the emphasis was on reducing avoidable emergency attendances or admissions, for conditions such as asthma, urinary tract infections and chronic obstructive pulmonary disease that are treatable in a primary setting given timely and appropriate care [2]. In the last few years, the UK's National Health Service (NHS) strategy has switched to containment of emergency care demand in general, with the government seeking to reduce unplanned emergency admissions to hospital for long-term chronic conditions [3,4], for example, by greater use of community matrons [5] assigned to patients with high intensity health needs and by encouraging primary and community care options for such patients wherever possible. Such initiatives are sometimes grouped as constituting the "chronic care model" [6]. Their cost effectiveness has been demonstrated for appropriately selected conditions [7] -that is, care is cheaper without any deterioration in clinically defined measures of patient health state [8].
GP practices have a major role in managing the demand generated by long term conditions but may vary in their referral rates for avoidable emergencies [9] or have access arrangements that leave patients with little choice but to resort to hospital emergency units [10,11]. Some studies have concluded that there are significant variations in practice referral behaviour [12,13], and the UK government has identified wide variations [4] in rates of unplanned emergency admissions between health areas. Variations in referrals unrelated to need may in part be related to resourcing mechanisms that do not sufficiently compensate for variations in ill health between population subgroups [14]; spatial access issues are also important [15]. Improved resourcing of, and access to, primary care in deprived areas means primary and community care teams will be better able to address patient morbidity in community settings, and patients will be less likely to attend hospitals for minor emergencies.
The present paper develops methods that recognises the interplay between population mix factors (e.g. social composition, health need and genuine morbidity differences), access, and primary care factors, that allows for unobserved influences on hospitalisation usage, and that also reflects interdependence between responses. Tools such as these may assist health care agencies to take actions to improve equity in access to primary care, to target chronic care model initiatives and to improve the match of resources to need, with the aim of reducing unplanned emergency admissions or attendances at hospitals.

Contextual Setting
In the UK health system, patients may choose their GP practice and there are no geographic constraints on their choice, so patients resident in small areas i = 1,...,I are typically affiliated to a range of practices (j = 1,..,J), though in practice geographic access to the GP surgery strongly affects choice of practice. This paper proposes a modelling strategy for multiple health outcomes (e.g. admissions, attendances) where "outcomes" is used in a broad statistical sense as a response variable rather than the health care sense. The strategy allows for the impact on referral rates of known factors such as area social structure (e.g. area deprivation), and differences in geographic accessibility to both primary care and to hospitals with emergency units. It also allows for unknown influences taken to vary randomly over areas and practices. This involves a cell decomposition approach whereby area populations (and the hospitalization or other health care events they generate) are disaggregated according to the primary care practices they are affiliated to. So variations for multiple outcomes k = 1,..,K are considered for population totals P ij classified both by area i and GP practice j.
Since in fact many cells in the I × J cross-classification of patients by area and practice are empty or contain very few patients (e.g. because practice j is too geographically remote from area i to attract any patients from it), the method of analysis in this paper uses a subset of the IJ possible area-practice intersections. This involves including area-practice cells which account for the great majority of a practice's population, but excluding those that may account for just a few patients. So the observations are totals Y kh by outcome k and cell h = 1,..,H, where H is less than IJ, and where i h and j h are the area and practice corresponding to cell h. Here K = 2, and attendances and admissions Y kh (k = 1 for attendances, k = 2 for admissions) are considered in relation to expected cases in the Various factors influencing attendance and admission rate variations between areas may be unobserved. The modelling strategy therefore includes an allowance for spatially correlated but unobserved risk factors for areas; these are modelled by an intrinsic conditional autoregressive or ICAR prior [16]. Similarly a practice level random effect summarises influences on practice referral rates that cannot be accounted for by known measures of practice health needs. Different approaches are possible to model correlated need and behaviour: for example, one may use a K-variate ICAR prior for area spatial effects combined with K-variate practice errors. Another option is a common factor approach which is less heavily parameterised (especially as K increases), and will be effective when there are high correlations in area and practice relativities over the outcomes.
The method extends the approach of Congdon and Best [17] to multivariate outcomes. It is adapted to a situation where multiple outcomes (e.g. attendances and admissions) are interrelated in terms of morbidity influences (variations in health need due to known and unknown factors) and also behavioural influences (e.g. referral thresholds of GP practices).
The methods described in the paper adopt a fully Bayesian strategy. This involves specification of prior densities on parameters and updating such densities (to provide posterior parameter densities) via the likelihood of the observed data. In estimating the models, iterative Monte Carlo Markov Chain (MCMC) techniques [18] are used, as implemented in the WINBUGS program [19]. The geographical unit of analysis is the recently introduced Super Output Area (abbreviated as SOA), designed by the UK Office of National Statistics Census Unit SOAs to avoid the problems caused by the inconsistent and unstable electoral ward geography and improve statistical comparison as they are of consistent size [21]. Specifically lower level SOAs are used here. There are I = 149 such SOAs in Havering with an average population P i of 1500. Figures 1 and 2 show the variation in maximum likelihood SOA attendance and admission ratios, namely standardised activity ratios (SARs) of observed to expected events. Figure 3 shows SOA indicators for deprivation rank (within Havering), as measured by Index of Multiple Deprivation (IMD) scores for 2004 [22]. These maps use quintile breaks.

Case Study
For the regression modelling of area and practice effects, the threshold approach mentioned above for including area-practice cells was applied with respect to SOA populations. If a GP practice accounted for more than 1% of an area's population the corresponding area-practice cell was included in the study. Under this definition, 97.7% of Havering's population who are affiliated to a Havering practice are included in the analysis, and there are H = 1620 cells.
Indirect standardisation is chosen in Figures 1 and 2 and in the regression analysis below because it is appropriate in the event of the small attendance/admission totals (including zeroes), and also often small population denominators. The regression analysis is at the area-practice population cell level and direct standardisation would involve cell level rates disaggregated by age and sex, obtained by dividing cell level events in each of 38 age-sex categories by cell level populations in each category. There SARs for Emergency Attendances Figure 1 SARs for Emergency Attendances.  are over 1600 cells in a population of around 220 thousand, so the population denominators involved are often small and the rates obtained have high sampling variability.

Standard Attendance Ratio
To quote the Clinical Indicators Support Team of NHS Scotland [23], "direct standardisation is inadvisable if the number of cases in any of the cells of the cross-classification of the variables used to standardise is small. Thus if one is standardising for age, sex and deprivation and there is a possibility of very low numbers in any combination of age, sex and deprivation categories, direct standardisation should be avoided. If there is a possibility that there are no cases in any of the cells of classification (zero cells) then direct standardisation is entirely ruled out. Indirect standardisation is highly robust in the context of small cell numbers". A further feature of indirect standardisation reported by Kendrick & McCloed [24] "is that..it is difficult to envisage a more robust technique for case mix adjustment. Not only does it adjust for the effects of individual variables, such as age, on outcome, it automatically adjusts for any interactions resulting from particular combinations of case mix variables. Indirect standardisation literally adjusts for the effects of all combinations of the casemix variables. ... indirect standardisation produces exactly the same results (is logically equivalent to) standardisation based on a fully saturated logistic regression model".
In the analysis below the implementation of indirect standardisation can be represented generically by considering the outcome Y s for stratum s (e.g. area, practice or cell) as a Poisson variable with mean µ s = E s ρ s , where E s is the expected outcome total or "exposure" for a Poisson response; for a fuller discussion of this method see, for example, [25] and [26]. Written out as a log link regression one has log(µ s ) = log(E s ) + log(ρ s ) and since the implicit coefficient on log(E s ) is 1 the E s term is sometimes called an offset [27].
In the present application the expected events E s are based on internal standardisation using data from the GP Population Register on the age-sex structure of the population at risk, and with the standard rates based on the entire Havering population included in the analysis. When expected outcomes are calculated using internal standardisation in this way, the ρ s can be interpreted as relative risks with average 1.

Patterns in Maximum Likelihood Activity Ratios
There are strong associations between standardised attendance and admission rates and social deprivation (IMD scores for 2004) in the case study population [22]. Such scores can also be calculated for GP practices by averaging over the scores for their patient places of residence. Table 1 shows the highest all ages attendance and admission ratios are in deprived SOAs and practices. There are positive correlations (of 0.65 and 0.72 respectively) between maximum likelihood SOA attendance and admission ratios (ratios of observed to expected events, Y ik /E ik ) and the IMD scores for such SOAs.
One may apply significance tests (from a classical rather than Bayesian perspective) to the practice or area SARs or functions of them to assess possible outliers. Suppose a high referring practice is defined as one with a SAR exceeding 150 and with the 2.5% confidence point exceeding 100. Using the procedure proposed in [28], five practices (11,19,35,42,49) have standard attendance ratios and 95% intervals on them fulfilling this criterion, and four practices (11,35,36,50) have 95% standard admission ratios fulfilling the criterion. Practices 11 and 35 are outliers for both attendances and admissions.
However, such fixed effects maximum likelihood estimates (and tests based on them) are often based on relatively small event numbers and also do not take account of spatial correlation in adjacent area rates or of the potential for pooling inferential strength by considering the underlying population density of rates over areas, practices, or cells. As Leyland and Davies [29] mention, "a map displaying [maximum likelihood fixed effect] SMRs tends to be dominated by areas with small populations since small changes to the observed number of deaths will result in large changes to the SMR". The approach of the present paper is based on the principle of pooling strength over related units (such as areas, schools, hospitals or GP practices) by modelling the density underlying the popu- lation of effects, after controlling for known influences on such effects. This approach is well known in applications in disease mapping and meta-analysis -for a fuller discussion, see for example [30] and [31].
Random effects may also be needed to model excess heterogeneity or overdispersion [32] relative to that expected under the Poisson model -this is measured by comparing the scaled deviance [33] with the number of observations. As noted in [34], "One of the most common reason for data being over-dispersed is that experimental conditions are not perfectly under control, and thus the unknown parameters vary not only with measured covariates but with latent and uncontrolled factors." Nevertheless a fixed effects model based on measured covariates only is included in the case study as a baseline model.

Models for Geographical and Practice Variations: Fixed Effects Model based on Known Influences
The analysis of this paper is predicated on looking at both area (SOA) and practice effects simultaneously via a cell level analysis. This is argued to be relevant for a crossed analysis of hospital outcomes and provides control for the catchment area effects on a GP practice's workload (e.g. a practice has a high referral rate because it catchment area has a population with high morbidity). Conversely, a cell analysis also controls for the possibility that area effects are distorted by possibly unusual practice referral rates (e.g. when a GP practice in a low morbidity area has an unusually high referral propensity, due possibly to varying referral thresholds). By comparison, analysis at practice level alone, or at SOA level alone is partial, and neglects area-practice interactions. (Analogous problems of ecological or atomistic distortions occur when a single rather than multilevel analysis is performed with nested as opposed to crossed data). Therefore an appropriate regression analysis is at the area-practice cell level.
The model used here for event k totals from cell h (admissions Y 1h and attendances Y 2h ) takes them to result from risks of the event in the area i h and in the GP practice j h that define cell h. Assuming Poisson variation in the counts where E kh are expected admissions/attendances based on applying the study area wide (i.e. Havering) age-sex rates to the cell populations and ρ kh is the relative event rate in cell h (with mean 1 if internal standardisation is used). Define where α k is a constant term.
One possible approach to analyzing variations in referral rates over cells is to use Poisson regression (because the responses are counts) but without allowing for unobserved area or practice influences. Thus the regression would use only known attributes of areas, practices or cells in a fixed effects model.
The first such class of indicators are needs scores for practices and areas. The models below are framed in terms of composite area and practice health need scores. There are several potential indicators of need at area and/or practice level such as the IMD2004 total deprivation score, and more specific domain scores: namely the income domain score (the proportion of the population experiencing income deprivation in an area), and the "health, deprivation and disability" (HDD) domain. Another possibility is the 2001 Census standard illness ratio [35]. Such scores are calculated for GP practices by averaging over the scores for their patient places of residence. However, correlations between these indicators are high within the SOAs and within practices (see Table 2). Including (say) as independent variables both the SOA income scores and SOA health domain scores together with the practice level income and health domain scores might seem to be representing subtly different aspects of need, but would be at the expense or distorting the regression results. For simplicity therefore, the model specification is framed in terms of a single needs score for SOAs, and a single needs score for practices. These are denoted AHN i and PHN j respectively and have mean 0 and variance 1. They are based on maximum likelihood factor analysis (in SPSS) of the four indices in Table 2. The eigenvalues from the correlation matrix show overwhelming support for there being only a single underlying need score for areas and for practices.
Other known influences on hospital utilisation at area level include the Euclidean distance to the case study hospital (R i ), and access to primary care A i . The access score A i is based on the Euclidean distances d ij between SOAs and the GP practices j = 1,..,J (J = 53) that are based in Havering. This score takes account of the number M j of GPs in each GP practice. Thus where f(d) is a declining function of the distance d ij between the population centroid of area i and the location of the surgery of practice j. Where a practice has a branch surgery (5 out of 53 practices) the distance calculation is to the branch surgery. Here exponential decay is assumed, in line for instance with the health resource analysis by [36], such that f(d ij ) = exp(-bd ij ) b>0 (Equation 4) It is expected that an analysis based on travel times as opposed to Euclidean distance would give similar results. Use of travel times is problematic as visitors to GP surgeries may use one of several alternative modes, and additionally Euclidean distances tend to be strongly correlated with travel times [37].
It is also necessary to take account of the impact on celllevel attendance and admission rates of variations in primary care referral behaviours. As well as the practice's overall population morbidity level represented in the composite need score PHN j , the proportion NH j of the elderly (over 75) practice population that are living in nursing and residential homes may be relevant. These homes contain high proportions of frail elderly and are subject to high hospitalisation rates.
Other known practice predictors that may influence referrals were included in exploratory work (whether the practice is single handed and list size to GP ratio) but were not significant. These are also not indicators of health need per se, but of practice organisation. The goal here is to control for the impact of health need on practice behaviour and assess residuals as indicative of variations in referral or attendance unrelated to need. It is important to contextualise GP practice variations for known population casemix factors as relatively good or bad 'performance' indicators for (say) health providers or schools may result in part from the social structure and health need of their population [38].
Then a fixed effects approach using known influences only will be based on SOA need, primary care access, and distance to hospital, and practice need and nursing home patients. So for cell h and outcome k the fixed effects regression is

Random Effects Options
Even after accounting for known impacts on need and access to care, there are likely to remain unknown factors influencing hospital usage outcomes. Consider unknown area and practice influences on hospitalisation rates. The contribution of the present paper lies particularly in specifying random effects for such unknown factors in analysis of attendances Y 1h and admissions Y 2h in cell h, that is in SOA i h and practice j h , while taking account of relevant known area and practice influences. While classical estimation of random effects models has progressed in recent years, the view taken in the current paper is that Bayesian estimation considerably facilitates the estimation of relatively complex random effects models [39].
Unknown influences may be correlated over areas, outcomes or both. For example, neighbouring areas may experience similar impacts on morbidity (and hence hospitalisation) from unknown environmental factors. These effects are also likely to be correlated across outcomes, implying the need for a multivariate spatial error. For instance, densely populated areas close to major traffic routes may show both high attendance and admission rates for conditions related to air and traffic pollution. For practices, unobserved variations in practice referral behaviour are likely to be correlated over outcomes.
Consider first a K-dimensional error structure for both area and practice random effects. Let s ki denote spatially correlated influences on morbidity for outcome k and areas i. Then under the MCAR prior, S i has a conditional prior density where N K denotes a K dimensional normal density, S [i] denotes spatial effects other than those in area i, ρ∈[0,1] is a scalar, Σ s is a K × K covariance matrix, and W ij = w ij I K × K is K × K with w ij denoting row standardized spatial interactions, namely . where the remaining parameters are assigned fixed effect priors. Estimates of the log relative risk of event k for area i may be obtained by monitoring the quantities η (a) ki centred round their average over all SOAs. Posterior means of the exponentials of centred η (a) ki can be obtained to provide estimates µ ki of relative risk in area i.
There will also be many unobserved influences on primary care effectiveness that are here initially summarised in a practice level random effect for each outcome e kj . Then a K-variate error model for the GP practice component of η kh includes the composite practice health need score and nursing home effects, and errors e kj . So with j = j h η (p) kj = β k1 PHN j + β k2 NH j + e kj (Equation 9) where e j = (e 1j , e 2j ) are bivariate normal random effects with mean zero and covariance Σ e . So with K = 2 (e 1j , e 2j ) ~ N 2 (0, Σ e ), j = 1,..., j. (Equation 10) Posterior means of the exponentials of centred η (p) kj can be obtained to provide estimates of relative risk ν kj in practice j.
If a model based on practice and area predictors and random effects leaves substantial unexplained variation (i.e. produces a model from which predictions do not match the observations), it may be necessary to add predictors and/or random effects at cell level. For example, multivariate cell level random effects u h = (u 1h , u 2h ,..u Kh ), h = 1,...H, may be needed to represent remaining heterogeneity, as in where (u 1h , u 2h ,..u Kh ) ~ N K (0, Σ u ). It may be that only one outcome (say k*) is not satisfactorily predicted and in this case univariate unstructured effects u k*h at cell level may be introduced for that outcome only. Another possibility for model extension is heterogeneity over SOAs in the effects of known predictors in equation (8), such as varying impacts of GP access or distance to hospital. Such heterogeneity can itself be modelled using an MCAR prior over outcomes k.
Multivariate correlation in errors does not necessarily require a conventional K dimension error density (e.g. a K-variate normal). One possible alternative approach assumes that the s ki and the e kj are generated by a single area factor s i and a single practice factor e j respectively, with outcome specific loadings governing the relative importance of the common factor to explaining each outcome. Under a single common factor model the area and practice components have the form Equation 13) where {κ k , λ k } are loadings. Whether all the event specific loadings κ k and λ k are free parameters depends on whether the variances and of s i and e j are preset [41].

Policy Relevant Model Outputs
Measures to reduce high attendance rates may involve provision of extra primary care clinics or staff and the siting of such resources is important. A possible index of poor access in relation to health care need is the discrepancy between the attendance or admission relative risks in different areas (µ ki ) and levels of access to primary care.
The latter is measured by the ratio of access in SOA i to average access in all SOAs, namely A i / . So is a measure suggesting where demand for A&E attendances or admissions might be reduced by improving access to primary care. From equations 3 and 4, A i is stochastic so the ∆ ki need to be monitored during the MCMC sampling and a posterior density obtained. High posterior means in ∆ ki (or in the ranks of ∆ ki ) in particular SOAs will occur when a high attendance or admission rate is combined with access below average.
Also relevant are measures of practice referral behaviour after controlling for known influences on health need (i.e. controlling for the level of morbidity of the population that a GP practice serves). These are represented here by the residuals e kj under the multivariate error models, or the pooled e j under the common factor model. Posterior densities of ranks on such residuals are useful in assessing performance in terms of referral rates in excess of health need, and have been used in other institutional comparisons [42,43].

Model Estimation and Priors
The analysis starts by comparing two models, the first based on fixed effects only using the observed predictors. Thus for the fixed effects model The second model assumes a full dimension error structure for both practice and SOA effects but has no cell level effects. Relatively diffuse N(0,1000) priors are used for the fixed effects {α k , β k1 , β k2 , γ k1 , γ k2 , γ k3 } under all models.
For the full dimension error models, Wishart priors with identity scale matrix and K = 2 degrees of freedom are assumed for the precision matrices and .
The remaining unknown in both models is the distance decay parameter in the exponential model for GP access scores in equation ( Fit is measured using the Deviance Information Criterion (abbreviated DIC) of Spiegelhalter et al [44], considered in terms of the separate outcomes. Denote the posterior mean of the deviance as k = -2 k where L k is the loglikelihood for event k taken over all cells h. The DIC for outcome k is k plus a complexity penalty c k , estimated as c k = k -D( k ), where D( k ) is the deviance at the mean of the parameter set θ k for model k. The extent of overdispersion is assessed by obtaining the posterior mean of the scaled deviance.
The ability of the model to reproduce the data is assessed via predictive checks from the posterior predictive density P(Y new |Y), where model predictions (replicate data sampled from the model) are sampled at iterations t = T 1 ,..T 2 from the MCMC sample chains, and iterations 1,..T 1 constitute the convergence stage. Following Gelfand [45], a check is made whether the observed Y kh are within 95% intervals of the predictions Y new,kh . For a satisfactory model the 95% intervals of the predictions should include the actual observations with probability of 0.95 or higher.
Estimation uses two chains with dispersed initial values run for 10000 iterations. Convergence was assessed by Gelman-Rubin criteria [46], and summaries of the parameters, the discrepancy indices ∆ ki , and GP practice error ranks, are based on iterations 5,000-10,000. Table 3 shows that a random errors interdependent outcomes model (model 2) has better fit than the fixed effects model 1 as measured by the DIC, albeit at the cost of an increase in complexity. The fixed effects model shows overdispersion in relation to the total number of observations (namely 1610 for each outcome and 3220 for both outcomes combined). The random effects model provides satisfactory predictions of the data for attendances, and removes overdispersion in attendances. However, some predictive discrepancies remain for admissions, since only 91.4% of the 95% credible intervals for replicate data include the actual data. The scaled deviance for admissions under model 2 still indicates overdispersion.

Results
Coefficient summaries (Table 4) for these two models show that allowing for unobserved area effects (e.g. unobserved influences on morbidity and access) enhances the distance effect to hospital parameters, γ k3 . Such enhancements of fixed effect parameters (in absolute terms) are often obtained in random effects models. The coefficients  Model 3 is therefore used to make inferences for particular areas or practices that are of interest in health strategy terms. For practices, discrepancies against anticipated utilisation are assessed by the residuals e kj in equation (9), which control for known need and access influences on utilisation. Table 5 lists posterior medians and 95% intervals for ranks on residual practice effects e kj under model 3 (highest positive residuals have higher ranks), together with practice population needs scores, obtained via factor analysis from four measured indicators. Also provided are estimated practice relative risks ν kj by event k.
There is a correlation of 0.74 between the posterior ranks on e 1j and those on e 2j , demonstrating overlapping issues of GP practice referral behaviour for the two outcomes. There is overlap between the results in Table 5 and the earlier reported analysis based on maximum likelihood SARs and classical 95% intervals. However, practices 36 and 50 no longer figure as leading admission rate outliers (in terms of median ranks on residuals e k2 ) after taking account of their need scores and nursing home patients; similarly, practices 19 and 49 no longer figure as attendance rate outliers.
Outlier areas are also important to identify when excess hospital attendances or admissions by small area coincide with relatively poor access to primary care. Figures 4 and  5 shows the distribution of the median posterior ranks of the discrepancies in equation (14). These two sets of discrepancies are highly related (a correlation of 0.94 between the 149 posterior medians for the two outcomes). There are also positive correlations between such discrepancies (as measured by posterior median ranks) and the event relative risks themselves (0.62 for attendances, and 0.59 for admissions). There are also positive correlations between the discrepancies in equation (14) and area IMD scores (0.49 for attendances, and 0.30 for admissions), demonstrating an interaction between relatively low primary care access and social deprivation in causing potentially avoidable hospital utilisation.

Discussion
Modern health care strategy emphasises the need for avoiding unnecessary referrals to acute care and the need for effective health care strategies for long term chronic conditions that are based in primary and community care settings. The paper has sought to present a methodology that analyses interactions between population casemix and the referral behaviour of primary care "gatekeepers" and that also reflects interdependence between outcomes. Such a methodology can be widely applied as ensuring an effective and appropriate balance between primary and secondary care in the management of chronic disease is a recurrent issue. Particular goals might include identifying whether there are areas or GP practices with significantly elevated referral rates, especially when not justified by morbidity.
The analysis of the paper has had a specialised focus on a relatively affluent outer city area and on two particular outcomes, namely emergency admissions and attendances at hospitals. Such forms of hospital use are often assessed as unnecessary, especially when involving conditions that could be treated in primary care settings. The present analysis identified GP practices that had high attendance and admission rates in relation to the needs level in their catchment populations (e.g. practices 11, 35 and 42 in Table 5), and this is also apparent in median ranks on residuals after controlling for need. While admitting that health care need may often have multiple separate dimensions, in the case study a selection of commonly used indicators seem to be measuring essentially the same dimension (Table 2) and a factor analysis was used to combine the multiple measured needs indices.
The present study has considered multiple interrelated outcomes but only for a single hospital. For many populations (e.g. in rural areas, relatively small towns, or outer suburban areas as in the case study) this reflects a reality of a single effective provider except for less common elective procedures where 'out of area' referrals are possible. However, many urbanised inner city areas have a choice of providers (e.g. general hospitals with emergency units).
If observations Y khm are available by area-practice cell h, outcome k, and relevant providers m = 1,..M, then the model can be provider specific, with Y khm ~ Po(E kh ρ khm ). This model would include supply factors (e.g. bed/staff numbers) in an additional provider/facility term η (f) km , as well as allowing for residence-hospital interactions η (a) kim . So for i = i h and j = j h In the origin-hospital interaction component, distances to hospital (R im ) would be specific to each area and provider, so that for i = i h η (a) kim = γ k1 AHN i + γ k2 A i + γ k3 R im + s ki . (Equation 19) Alternatively, an access to hospital term (allowing for both hospital mass and distances from areas to hospitals) could replace R im in equation 19.
The present study has considered particular relatively broad outcomes and a particular geographic setting. Different outcomes, specific to particular diagnostic or disease groups, or to different types of geographic setting, would need to adapt the model. For instance, compared to other London boroughs, Havering is relatively affluent, with lower illness and deprivation levels than average. However, it is internally diverse, containing some deprived pockets of social rented housing. In other London boroughs features such as ethnic mix may play a greater role in explaining some types of referral activity [47]. Havering is a largely affluent and white borough (around 95% white in 2001, the highest in London), with some resemblance to less metropolitan areas outside London, in having pockets of deprivation, and in having low homeless and refugee populations. Analysis of particular outcomes in different geographic settings, such as cardiovascular or mental illness referrals in inner city areas, would need to include ethnicity, as well as the significant presence of marginalised groups, as influences on need.