Skip to main content

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



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.


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.


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.


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 Pij 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 Ykh by outcome k and cell h = 1,..,H, where H is less than IJ, and where ih and jh are the area and practice corresponding to cell h. Here K = 2, and attendances and admissions Ykh (k = 1 for attendances, k = 2 for admissions) are considered in relation to expected cases in the population living in area ih and affiliated to practice jh.

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].

Case Study

A case study analysis uses data on 20200 attendances and 5970 emergency admissions at Oldchurch Hospital in North East London during April 1 – September 30th 2003, and focuses on attendances at, and admissions to, this hospital by residents of the outer London borough of Havering, within which the hospital is located. The majority (93.7%) of the 239 thousand residents of Havering (according to the NHS patent register) are affiliated to one of J = 53 GP practices sited within Havering; around 15000 are affiliated to practices outside Havering.

The data sources are standard NHS patient recording systems. The attendances data are based on a patient administration system maintained by the Barking, Havering and Redbridge Trust (BHRT), while the admissions data are from the centrally collated NHS Health Episode Statistics, widely used (e.g. in monitoring government health targets) to analyse morbidity, access and utilisation in the English health care system. Oldchurch accounted for around 80% (5970/7510) of emergency admissions by Havering residents in the period. It is not possible to ascertain how many attendances there are by Havering residents to all hospitals in the surrounding region, since attendance data are not collated centrally, and individual trusts have their own systems (or may not have attendance recording systems at all); attendances by Havering residents at Oldchurch greatly exceed those at King George, the other BHRT hospital with an emergency unit. Patient populations by age and sex are obtained from the National Health Service Central Register [20] which holds demographic details for patients registered with General Practitioners in different health areas – currently such areas are known as Primary Care Trusts (PCTs) and Havering PCT is coterminous with the London Borough of Havering. Census data on long term illness are from the 2001 Census, as discussed below.

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 Pi 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.

Figure 1
figure 1

SARs for Emergency Attendances.

Figure 2
figure 2

SARs for Emergency Admissions.

Figure 3
figure 3

IMD2004 Multiple Deprivation Score.

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 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.

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 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 Ys for stratum s (e.g. area, practice or cell) as a Poisson variable with mean μs = Esρs, where Es 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(Es) + log(ρs) and since the implicit coefficient on log(Es) is 1 the Es term is sometimes called an offset [27].

In the present application the expected events Es 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, Yik/Eik) and the IMD scores for such SOAs.

Table 1 Maximum Likelihood Activity Ratios (Admissions & Attendances) by Deprivation Category, SOAs and GP Practices

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 population 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 Y1h and attendances Y2h) takes them to result from risks of the event in the area ih and in the GP practice jh that define cell h. Assuming Poisson variation in the counts

Ykh ~ Po(Ekh ρkh),     k = 1, ... K     (Equation 1)

where Ekh 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

log(ρkh) = αk + ηkh     (Equation 2)

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.

Table 2 Correlations of needs scores* over practices & SOAs

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 AHNi and PHNj 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 (Ri), and access to primary care Ai. The access score Ai is based on the Euclidean distances dij between SOAs and the GP practices j = 1,..,J (J = 53) that are based in Havering. This score takes account of the number Mj of GPs in each GP practice. Thus

where f(d) is a declining function of the distance dij 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(dij) = exp(-bdij)     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 cell-level 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 PHNj, the proportion NHj 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

ηkh = γk1AHNi + γk2Ai + γk3Ri + βk1PHNj + βk2NHj.     (Equation 5)

This model can be estimated by classical techniques (e.g. in STATA or R), or by Bayesian techniques. Except for small samples, classical and Bayesian estimates of such a model are usually very similar.

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 Y1h and admissions Y2h in cell h, that is in SOA ih and practice jh, 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 ski denote spatially correlated influences on morbidity for outcome k and areas i. For modelling correlation between ski and smi (k ≠ m), Gelfand and Vounatsou [40] propose a multivariate conditional autoregressive (MCAR) prior. Let Si = (s1i,...sKi) be the multidimensional error for the ith area. Then under the MCAR prior, Si has a conditional prior density

where NK 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 Wij = wijIK × K is K × K with wij denoting row standardized spatial interactions, namely . If the wij are based on contiguity (that is, cij = 1 if SOAs i and j are adjacent, and cij = 0 otherwise), then wij = 1/Li if areas i and j are adjacent, with Li being the number of neighbours of area i. The introduction of ρ ensures the corresponding joint prior is proper, with nonsingular covariance matrix. If ρ is set to 1, as in

then a propriety issue occurs, though identifiability may be achieved by centering each of the K sets of effects at each MCMC iteration. This approach is used in WINBUGS and in the analyses reported here.

Then the model for the area component of ηkh involves area ih that together with practice jh defines cell h. So for i = ih

η(a)ki = γk1AHNi + γk2Ai + γk3Ri + ski     k = 1,...,K     (Equation 8)

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 ekj. 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 ekj. So with j = jh

η(p)kj = βk1PHNj + βk2NHj + ekj     (Equation 9)

where ej = (e1j, e2j) are bivariate normal random effects with mean zero and covariance Σe. So with K = 2

(e1j, e2j) ~ N2(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 uh = (u1h, u2h,..uKh), h = 1,...H, may be needed to represent remaining heterogeneity, as in

where (u1h, u2h,..uKh) ~ NK(0, Σu). It may be that only one outcome (say k*) is not satisfactorily predicted and in this case univariate unstructured effects uk*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 ski and the ekj are generated by a single area factor si and a single practice factor ej 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

η(a)ki = γk1AHNi + γk2Ai + γk3Ri + κksi     (Equation 12)

η(p)kj = βk1PHNj + βk2NHj + λkej     (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 si and ej 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 Ai/. So

Δki = μki-Ai/     (Equation 14)

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, Ai 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 ekj under the multivariate error models, or the pooled ej 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

η(a)ki = γk1AHNi + γk2Ai + γk3Ri     (Equation 15)

η(p)kj = βk1PHNj + βk2NHj.

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 (4. To model b as a continuous unknown would be at the cost of considerably greater computing times, because the need for repeated calculation of f(dij) involves a 149 × 53 distance matrix. So a 24 point discrete prior for b on values {0.1,0.15,0.2,0.25,0.3,...,0.9,0.95,1,1.1,1.2,1.3,1.4,1.5} is assumed using pre-calculated decay matrices at these values. This range of values accords with prior knowledge; for example, Carr-Hill et al [36] assume b = 0.2.

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 = -2k where Lk is the log-likelihood for event k taken over all cells h. The DIC for outcome k is k plus a complexity penalty ck, estimated as ck = 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(Ynew|Y), where model predictions (replicate data sampled from the model)

are sampled at iterations t = T1,..T2 from the MCMC sample chains, and iterations 1,..T1 constitute the convergence stage. Following Gelfand [45], a check is made whether the observed Ykh are within 95% intervals of the predictions Ynew,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.

Table 3 Model Fit and Checking Criteria

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 γk1 and βk1 reflecting the impact of composite health need are both significant, and enhanced in model 2. The impacts of access to primary care (increases in which might be expected to reduce emergency admissions and attendances), as summarised in the parameter γk2, are negative. In both models the posterior densities for bk in fk(dij) concentrate on small values.

Table 4 Parameter Estimates

A third model seeks to improve fit and predictions for admissions by adding a cell level random error for k = 2. So

where u2h ~ N(0, ), and 1/ is assigned a gamma G(1,0.001) prior. Predictive checks for this model are satisfactory, and despite increased complexity for the admissions component (with c2 = 596), the DIC for all parts of the model is improved over model 2. The remaining overdispersion in admissions is only slight.

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 ekj 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 ekj 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.

Table 5 Practice Level Outcomes from Area-Practice Model 3

There is a correlation of 0.74 between the posterior ranks on e1j and those on e2j, 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 ek2) 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.

Figure 4
figure 4

Ranks on Attendance Residuals.

Figure 5
figure 5

Ranks on Admission Residuals.


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 Ykhm are available by area-practice cell h, outcome k, and relevant providers m = 1,..M, then the model can be provider specific, with Ykhm ~ Po(Ekhρ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 = ih and j = jh

In the origin-hospital interaction component, distances to hospital (Rim) would be specific to each area and provider, so that for i = ih

η(a)kim = γk1AHNi + γk2Ai + γk3Rim + ski.     (Equation 19)

Alternatively, an access to hospital term (allowing for both hospital mass and distances from areas to hospitals) could replace Rim 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.


  1. British Columbia Medical Association: Ensuring Excellence: Renewing BC's Primary Care System. BCMA Policy Papers. 2002

    Google Scholar 

  2. Weissman J, Gatsonis C, Epstein A: Rates of avoidable hospitalization by insurance status in Massachusetts and Maryland. Journal of the American Medical Association. 1992, 268: 2388-2394. 10.1001/jama.268.17.2388.

    Article  PubMed  CAS  Google Scholar 

  3. Chaplin S: Cutting costs by reducing emergency admissions. Prescriber. 2006, 28-30.

    Google Scholar 

  4. Department of Health Press Notice: Improve healthcare by reducing unnecessary emergency admissions – Hewitt. Published: Monday 20 March 2006, Ref: 2006/0104

  5. Department of Health Press Notice: Community Matrons will make a difference to patients' lives – Chief Nursing Officer. Tuesday 1 February 2005. Ref: 2005/0036

  6. Singh D, Ham C: Improving Care for People with Long Term Conditions. University of Birmingham and NHS Institute for Innovation and Improvement. 2005

    Google Scholar 

  7. MacIntyre C, Ruth D, Ansari Z: Hospital in the home is cost saving for appropriately selected patients: a comparison with in-hospital care. Int J Qual Health Care. 2002, 14: 285-293. 10.1093/intqhc/14.4.285.

    Article  PubMed  Google Scholar 

  8. Donaldson C, Currie G, Mitton C: Cost effectiveness analysis in health care: contraindications. Brit Med Jour. 2002, 325: 891-894. 10.1136/bmj.325.7369.891.

    Article  Google Scholar 

  9. Saxena S, George J, Barber J, Fitzpatrick J, Majeed A: Association of population and practice factors with potentially avoidable admission rates for chronic diseases in London: cross sectional analysis. Journal of the Royal Society of Medicine. 2006, 99: 81-89. 10.1258/jrsm.99.2.81.

    Article  PubMed  PubMed Central  Google Scholar 

  10. Vedsted P, Christensen M: The effect of an out-of-hours reform on attendance at casualty wards: the Danish example. Scand J Prim Health Care. 2001, 19: 95-98. 10.1080/028134301750235303.

    Article  PubMed  CAS  Google Scholar 

  11. Carlisle R, Groom L, Avery A, Boot D, Earwicker S: Relation of out of hours activity by general practice and accident and emergency services with deprivation in Nottingham: longitudinal survey. Brit Med Jour. 1998, 316: 520-523.

    Article  CAS  Google Scholar 

  12. O'Donnell C: Variation in GP referral rates: what can we learn from the literature?. Family Practice. 2000, 17: 462-471. 10.1093/fampra/17.6.462.

    Article  PubMed  Google Scholar 

  13. Wilkin D, Smith A: Explaining variation in general practitioner referrals to hospital. Fam Pract. 1987, 4: 160-169.

    Article  PubMed  CAS  Google Scholar 

  14. Gulliford M: Availability of primary care doctors and population health in England: is there an association?. J Public Health Medicine. 2002, 24: 252-254. 10.1093/pubmed/24.4.252.

    Article  Google Scholar 

  15. Guagliardo M: Spatial accessibility of primary care: concepts, methods and challenges. International Journal of Health Geographics. 2004, 3: 1-13. 10.1186/1476-072X-3-3.

    Article  Google Scholar 

  16. Besag J, York J, Mollie A: Bayesian image restoration with two applications in spatial statistics. Annals of the Institute of Statistical Mathematics. 1991, 43: 1-59. 10.1007/BF00116466.

    Article  Google Scholar 

  17. Congdon P, Best N: Small area variation in hospital admission rates: adjusting for referral and provider variation. J Roy Stat Soc, Series C. 2000, 49: 207-226. 10.1111/1467-9876.00188.

    Article  Google Scholar 

  18. Gelfand A, Smith A: Sampling-based approaches to calculating marginal densities. J American Statistical Association. 1990, 85: 398-409. 10.2307/2289776.

    Article  Google Scholar 

  19. Spiegelhalter D, Thomas A, Best N, Lunn D: WinBUGS User Manual, Version 1.4. 2003, []

    Google Scholar 

  20. Hedley A, McMaster W: Use of the National Health Service Central Register for medical research purposes. Health Bulletin. 1988, 46: 63-68.

    PubMed  CAS  Google Scholar 

  21. ONS Neighbourhood Statistics. []

  22. Noble M, Wright G, Dibben C, Smith G, Maclennan D, Anttila C, Barnes H: The English Indices of Deprivation. Office of the Deputy Prime Minister, Neighbourhood Renewal Unit. 2004

    Google Scholar 

  23. Clinical Indicators Support Team, NHS Scotland. []

  24. Kendrick S, MacLeod M: Adjusting outcomes for case mix: indirect standardisation and logistic regression. CIST Working Paper No 3. []

  25. Wakefield J, Best N, Waller L: Bayesian approaches to disease mapping. Spatial Epidemiology: Methods and Applications. Edited by: Elliott P, Wakefield J, Best N, Briggs D. 2000, Oxford University Press

    Google Scholar 

  26. Clayton D, Kaldor J: Empirical Bayes estimates of age-standardized relative risks for use in disease mapping. Biometrics. 1987, 43: 671-81. 10.2307/2532003.

    Article  PubMed  CAS  Google Scholar 

  27. Wikipedia on Poisson Regression. []

  28. Ulm K: A simple method to calculate the confidence interval of a standardized mortality ratio. American Journal of Epidemiology. 1990, 131 (2): 373-375.

    PubMed  CAS  Google Scholar 

  29. Leyland A, Davies C: Empirical Bayes methods for disease mapping. Statistical Methods in Medical Research. 2005, 14: 17-34. 10.1191/0962280205sm387oa.

    Article  PubMed  Google Scholar 

  30. Spiegelhalter D, Abrams K, Myles J: Bayesian Approaches to Clinical Trials and Health Care Evaluation. 2004, Wiley, New York

    Google Scholar 

  31. Berke O: Exploratory disease mapping: kriging the spatial risk function from regional count data. Int J Health Geogr. 2004, 3: 18-10.1186/1476-072X-3-18.

    Article  PubMed  PubMed Central  Google Scholar 

  32. Dean C: Testing for overdispersion in Poisson and binomial regression models. Jour Amer Stat Assoc. 1992, 87: 451-457. 10.2307/2290276.

    Article  Google Scholar 

  33. McCullagh P, Nelder J: General Linear Models. 1989, Chapman & Hall, London

    Book  Google Scholar 

  34. Pedan A: Analysis of Count Data Using the SAS System. []

  35. Gordon D: Census based deprivation indices: their weighting and validation. Journal of Epidemiology and Community Health. 1995, 49: S39-S44.

    Article  PubMed  PubMed Central  Google Scholar 

  36. Carr-Hill R, Hardman G, Martin S, Peacock S, Sheldon TA, Smith P: A formula for distributing NHS revenues based on small area use of hospital beds. Centre for Health Economics Occasional Paper. 1994, York: University of York

    Google Scholar 

  37. Phibbs C, Luft H: Correlation of travel time on roads versus straight line distance. Medical Care Research and Review. 1995, 52: 532-542.

    Article  PubMed  CAS  Google Scholar 

  38. Goldstein H, Spiegelhalter D: League tables and their limitations: statistical issues in comparisons of institutional performance. J Royal Statistical Soc Series A. 1996, 159: 385-409.

    Article  Google Scholar 

  39. Best N, Spiegelhalter D, Thomas A, Brayne C: Bayesian Analysis of Realistically Complex Models. Journal of the Royal Statistical Society, Series A. 1996, 159 (2): 323-342.

    Article  Google Scholar 

  40. Gelfand A, Vounatsou P: Proper Multivariate Conditional Autoregressive Models for Spatial Data Analysis. Biostatistics. 2003, 4: 11-25. 10.1093/biostatistics/4.1.11.

    Article  PubMed  Google Scholar 

  41. Skrondal A, Rabe-Hesketh S: Generalized Latent Variable Modeling: Multilevel, Longitudinal and Structural Equation Models. 2004, Boca Raton, FL: Chapman & Hall/CRC

    Book  Google Scholar 

  42. Laird N, Louis T: Bayes and empirical Bayes ranking methods. Journal of Educational Statistics. 1989, 14: 29-46. 10.2307/1164724.

    Article  Google Scholar 

  43. C Guarino C, Ridgeway G, Chun M, Buddin R: A Bayesian latent variable model for institutional ranking. Higher Education in Europe. 2005, 30: 147-165. 10.1080/03797720500260033.

    Article  Google Scholar 

  44. Spiegelhalter D, Best N, Carlin B, van der Linde A: Bayesian measures of model complexity and fit. J Royal Stat Soc B. 2002, 64: 583-639. 10.1111/1467-9868.00353.

    Article  Google Scholar 

  45. Gelfand A: Model determination using sampling based methods. Markov Chain Monte Carlo in Practice. Edited by: Gilks W, Richardson S, Spieglehalter D. 1996, Chapman and Hall/CRC, Boca Raton, 145-161.

    Google Scholar 

  46. Gelman A, Carlin J, Stern H, Rubin D: Bayesian Data Analysis. 1995, Chapman and Hall, London

    Google Scholar 

  47. Morris S, Sutton M, Gravelle H: Inequity and inequality in the use of health care in England: an empirical investigation. 2003, Technical Paper Series 27, Centre For Health Economics, University of York

    Google Scholar 

Download references

Author information

Authors and Affiliations


Corresponding author

Correspondence to Peter Congdon.

Authors’ original submitted files for images

Rights and permissions

This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Reprints and permissions

About this article

Cite this article

Congdon, P. Modelling multiple hospital outcomes: the impact of small area and primary care practice variation. Int J Health Geogr 5, 50 (2006).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: