here - The BIAS project

Transcript here - The BIAS project

Models for small area data with applications in health care

Nicky Best

Department of Epidemiology and Biostatistics School of Public Health, Imperial College London http://www.bias-project.org.uk

http://www1.imperial.ac.uk/medicine/people/n.best/

Outline

1. Introduction 2. Mapping and spatial smoothing of health data 3. Classifying areas and summarising geographical variations in health outcomes 4. Modelling and mapping multiple outcomes 5. Modelling and mapping temporal trends 6. Hierarchical related regression models for combining individual and small area data

Outline

A brief history of disease mapping

• Health indicator maps have a long history in epidemiology and public health Spot maps: – Yellow fever pandemic New York (Seaman, 1798) – Cholera and the Broad Street Pump (Snow, 1854)

Spot map of cholera cases (Snow, 1854)

A brief history of disease mapping

Health indicator maps have a long history in epidemiology and public health • Spot maps: – Yellow fever pandemic New York (Seaman, 1798) – Cholera and the Broad Street Pump (Snow, 1854) • Chloropeth maps: – Geographical distribution of mortality from heart disease, cancer and TB in England & Wales (Haviland, 1878) – Cancer mortality by county in England & Wales, adjusted for age and sex (Stocks, 1936, 1937, 1939)

Female cancer 1851-60 (Haviland 1878)

Female lung cancer SMR 1921-30 (Stocks, 1939)

A brief history of disease mapping

Health indicator maps have a long history in epidemiology and public health • Spot maps: – Yellow fever pandemic New York (Seaman, 1798) – Cholera and the Broad Street Pump (Snow, 1854) • Chloropeth maps: – Geographical distribution of heart disease, cancer and TB in England & Wales (Haviland, 1878) – Cancer rates by county in England & Wales, adjusted for age and sex (Stocks, 1936, 1937, 1939) • National and international disease atlases, e.g

– Atlas of Cancer Incidence in England & Wales 1968-85 (Swerdlow & dos Santos Silva, 1993) – Atlas of Mortality in Europe 1980/81 & 1990/91 (WHO, 1997)

Female lung cancer incidence 1968-85 (Swerdlow and dos Santos Silva, 1993)

Age-standardised mortality from IHD, 1980-81 (WHO)

Recent developments in disease mapping

• Development of Geographical Information Systems (GIS) – Geographically indexed relational database – Computer program to map and analyze spatial data • Increasing availability of geo-referenced data – Ability to geocode, use GPS – Disease outcomes, demographics, environmental quality, health services • Development of statistical methods – Sophisticated techniques for separating signal from noise – Ability to account for spatial (and temporal) dependence – Methods for cluster detection and classification of areas Interest in mapping health events at small-area scale

Small area health data in the UK

• Administrative geography in UK includes – Postcodes (10-15 households) – Census Output Areas (COA; ~300 people) – Electoral wards (~500 to 2000 people) – Local authority districts, Health authority districts (10’s of thousands) • Postcoded data on mortality, births/still births, congenital anomalies, cancer incidence, hospital admissions • Population and socio-economic indicators from Census (COA) • Increasing availability of modelled environmental data at fine geographical resolution (grids) • Limited access to geographical identifiers for certain individual-level cohorts (e.g. Millenium Cohort, British Household Panel Survey) and health surveys (e.g. Health Survey for England)

Small area health data in Spain

• Administrative geography in Spain is divided into: – 17 regions – 52 provinces – ~8000 municipalities, ranging from small villages to large cities – Census tracts (finer sub-division in large cities) • Geocoded (place of residence) data on births, mortality (national), cancer incidence (regional; ~26% population), hospital discharge administrative data (national; public hospitals) • Small area (municipality) data on population and socio economic indicators from Census

Examples of recent disease atlases and health-related maps for Spain

Atlas of cancer mortality and other causes of death in Spain 1978-1992 ( López-Abente et al., 1996) Maps Age adjusted Rates and Standardised Mortality Ratios (SMR) at province level

Atlas of cancer mortality at municipality level in Spain 1989-1998 ( López-Abente et al., 2007) Maps (Bayesian) smoothed relative risks of mortality and probability of excess risk, at municipality level

Also produced maps of mortality from selected causes other than cancer, e.g. Influenza

…… + contextual maps of socioeconomic variables and environmental hazards

Outline

Why map small area disease rates?

• Interest in mapping geographical variations in health outcomes a the small area scale – Highlight sources of heterogeneity and spatial patterns – Suggest public health determinants or aetiological clues • Small scale – UK: electoral ward or census output area – SPAIN: municipality – – less susceptible to ecological (aggregation) bias better able to detect highly localised effects

Why smooth small area disease rates?

• Typically dealing with rare events in small areas A i Y i is the observed count of disease in area A i E i is the expected count based on population size, adjusted for age, sex, other strata …., • Relative risk usually estimated by SMR i = Y i / E i • Standard practice is to map SMRs

BUT sparse data need more sophisticated statistical analysis techniques

Why smooth small area disease rates?

• SMR represents estimate of ‘true’ (underlying) risk in an area, R i , i.e. R i = SMR i • Statistical uncertainty about estimate based on assuming Poisson sampling variation for data

Y

~ Poisson(R

E

) • SE( R

) = SE(SMR

)



1 / E

i – SMR i very populations imprecise for rare diseases and small – precision can vary widely between areas

Why smooth small area disease rates?

• SMR i in each area is estimated independently – ignores possible spatial correlation between disease risk in nearby areas due to possible dependence on spatially varying risk factors – leads to problems of multiple significance testing

Map of SMR of adult leukaemia in West Midlands Region, England 1974-86 (Olsen, Martuzzi and Elliott,

BMJ

1996;313:863-866).

Is the variability real or simply reflecting unequal expected counts ?

Have the red highlighted areas truly got a raised relative risk?

Methods for smoothing disease maps

• These problems may be addressed by spatial smoothing of the raw data • Idea is to “borrow information” for neighbouring areas to produce better (more stable, less noisy) estimates of the risk in each area • Similar principle to scatter plot smoothers, moving average smoothers….

• Many methods available

Methods for smoothing disease maps

• Ad hoc, local smoothing algorithms e.g. spatial moving averages, headbanging algorithm + quick and simple to implement – – can be very sensitive to ad hoc choice of weights etc.

no uncertainty estimates (standard errors) • Trend surface analysis e.g. kriging, polynomial/spline smoothing + estimation of ‘smoothing parameters’ based on trade off between fit and smoothness – + can be sensitive to choice of penalty for trade-off standard errors usually available

Methods for smoothing disease maps

• Random effects models e.g. empirical Bayes, hierarchical Bayes + data-based estimation of model parameters that + control smoothing full power of statistical modelling available: standard errors, prediction, probability calculations, inclusion of – covariates more complex to understand and implement

Bayesian Approach

• Use probability model to obtain smoothed risk estimate R i average) of in area i that is a compromise (weighted – observed area-level risk ratio ( Y i / E i ) – local or regional mean relative risk ( m ) • Weights depend on the precision (  1 / E i ) in area i and the of the SMR variability (heterogeneity) of the true risks across areas ( v ) local or regional mean relative risk ( m ) • Aim is to estimate posterior probability distribution of the unknown model parameters (R i , m, v) conditional on the data (Y i / E i )

Bayesian disease mapping model

• Typical Bayesian disease mapping model: Y i ~ Poisson(R i E i ), log (R i ) ~ Normal (m, v) • Hierarchical Bayesian model also requires specification of a (prior) probability distribution for m and v – These are often taken to be ‘non-informative’ • Empirical Bayes involves 2-step process: – Estimate m and v empirically from observed data – Ignore uncertainty in estimates of m and v and plug these values into the Bayesian model above

Software

• Estimation of Bayesian hierarchical models requires computationally intensive simulation methods (MCMC) – Implemented in free WinBUGS and GeoBUGS software: www.mrc-bsu.cam.ac.uk/bugs • Free software INLA (Rue et al, 2008) implements fast approximation: www.r-inla.org

• Empirical Bayes smoothing implemented in Rapid Inquiry Facility (RIF) : www.sahsu.org/sahsu_studies.php#RIF

Map of occurrences of adult leukaemia in West Midlands Region, England 1974-86 (A) unsmoothed SMR (B) smoothed by Bayesian methods (Olsen, Martuzzi and Elliott,

BMJ

1996;313:863-866)

Comparison of estimation methods

SMR Hierarchical Bayes Empirical Bayes 5 10 area 15 20

Including spatial dependence in disease risk

• R i are typically spatially correlated because they reflect, in part, spatially varying risk factors Incorporation of spatial dependence distribution of the R i ’s in the Conditional Autoregressive (CAR) model log (R i ) ~ Normal (m i , v i ) m i =  k R k / n i = average risk in neighbouring areas v i = v / n i → variance inversely proportional to number of neighbours Besag, York, Mollie (1991)

Annals of the Institute of Statistics and Mathematics

, 43: 1-59

RR N SMR N Raw data (SMR) Non-spatial smoothing (posterior mean R N i ) (418) < 0.5

(7) 0.5 - 0.7

<0.5

(55) 0.9 - 1.1

(39) 0.5 - 0.7

RR N 0.9-1.1

20.0km

Spatial smoothing (posterior mean R i ) (0) < 0.5

(26) 0.5 - 0.7

(193) 0.7 - 0.9

(292) 0.9 - 1.1

(252) 1.1 - 1.4

(90) 1.4 - 2.0

(20) >= 2.0

20.0km

(0) < 0.5

(39) 0.5 - 0.7

(192) 0.7 - 0.9

(264) 0.9 - 1.1

(251) 1.1 - 1.4

(96) 1.4 - 2.0

(31) >= 2.0

20.0km

Mapping uncertainty

• Mapping the mean posterior value of R i does not make full use of the posterior distribution • Map posterior SD • Map Probability ( R i > 1) Note – this is not the same as a classical p-value 1.0

1.5

Relativ e Risk Relative Risk, R i 2.0

2.5

RR N

Posterior SD of relative risk estimates

RR N Posterior mean relative risk N Posterior sd of relative risk (39) 0.5 - 0.7

(192) 0.7 - 0.9

(264) 0.9 - 1.1

(251) 1.1 - 1.4

(96) 1.4 - 2.0

(31) >= 2.0

N 20.0km

RR <0.5

20.0km

SD (0) < 0 .2

<0.2

(0) < 0 .2

(32 4) 0.2 - 0.4

(41 8) 0.4 - 0.6

(91 ) 0 .6 - 0 .8

(20 ) 0 .8 - 1 .0

(5) 1.0 - 1.2

(15 ) >= 1 .2

20.0km

RR N

Posterior probability that relative risk > 1

RR N Posterior mean relative risk N Posterior probability that relative risk > 1 (39) 0.5 - 0.7

(192) 0.7 - 0.9

(264) 0.9 - 1.1

(251) 1.1 - 1.4

(96) 1.4 - 2.0

(31) >= 2.0

RR values for pp 0.5-0.7

20.0km

Prob (16 3) < 0 .25

<0.25

(32 1) 0.2 5 - 0 .5

0.25-0.50

(16 3) < 0 .25

(32 1) 0.2 5 - 0 .5

(25 0) 0.5 - 0.7 5 (13 9) >= 0.7 5 20.0km

20.0km

Atlas of cancer mortality at municipality level in Spain 1989-1998 ( López-Abente et al., 2007)

Outline

Classifying areas with excess risk

• Richardson et al (2004): Simulation study investigating use of posterior probabilities in disease mapping studies • Classify an area as having an elevated risk if [Prob ( R i > 1)] > 0.8

• • High specificity (false detection < 10%) Sensitivity 60%-95% for E i 1.5-3.0 of 5-20 and true R i of

RR N

Childhood leukaemia in London

RR N Posterior mean RR N Posterior prob(RR>1) (39) 0.5 - 0.7

(192) 0.7 - 0.9

(264) 0.9 - 1.1

(251) 1.1 - 1.4

(96) 1.4 - 2.0

(31) >= 2.0

0.7-0.9

probability of RR greater than 1.0

20.0km

0.9-1.1

N <0.8

(77 2) < 0 .8

(10 1) >= 0.8

20.0km

Comparison of SaTScan and Bayesian classification rule SaTScan (Kuldorff, www.satscan.org): Bayesian: Location of most likely cluster N (772) < 0.8

(101) >= 0.8

< 0.8

2 nd 20.0km

Summarising geographic variation

• Often interested in providing overall summary measure of variability between areas, e.g.

– to compare variability of different outcomes – to quantify how much variation can be explained by covariates • Percentile Ratio: Ratio of outcomes (relative risks) in areas ranked at the q th and (100-q) th percentiles – e.g. 90 th Percentile Ratio, PR90 = R 95% /R 5% – Posterior distribution of PR90 easily calculated from MCMC output

Relative survival from colon cancer, England

Data • Survival/censoring times for all 7007 cases of colon cancer diagnosed in England in 1995 and followed for 5 years (provided by B Rachet, LSHTM) • Covariates: sex, age at diagnosis, clinical stage, deprivation score, Health Authority (95 area, 1-300 cases per HA) • Population mortality rates by age and sex for England and Wales, 1985 1995.

Questions of interest • Is there evidence of differences between Health Authorities in relative survival that may indicate differences in effectiveness of care received?

– Relative survival measures difference between age/sex-adjusted mortality rate in general population and in patients with disease of interest • How do these geographical differences change when we adjust for socioeconomic deprivation and clinical stage of cancer?

Relative survival from colon cancer, England

y kit ~ Poisson(m kit ) (subject k, area i, time interval t) log(m kit – E kit ) = log n kit + a t + bx ki + H i Area spatial effect Standard model for relative survival

Relative survival from colon cancer, England

y kit ~ Poisson(m kit ) (subject k, area i, time interval t) log(m kit – E kit ) = log n kit + a t + bx ki + H i Area spatial effect Standard model for relative survival N for deprivation and clinical stage relative excess hazard (model 6) N After adjustment for N Relative excess (9) 1.05 - 1.15

mortality (20) >= 1.15

<0.85

(7) < 0.85

(16) 0.85 - 0.95

(36) 0.95 - 1.05

(25) 1.05 - 1.15

(11) >= 1.15

PR90 = 1.95 (95 % CI 1.62-2.38) 200.0km

PR90 = 1.83 (95 % CI 1.54-2.24) 200.0km

200.0km

Ranking and classifying extreme areas

• Interest in ranking areas for e.g. policy evaluation, ‘performance’ monitoring • Rank of a point estimate is highly unreliable – Would like to measure uncertainty about rank • Straightforward to calculate posterior distribution of ranks (or any function of parameters) using MCMC – Obtain interval estimates for ranks • Can also calculate posterior probability that each area is ranked above a particular percentile

Rank (posterior mean and 95% CI) of the 95 Health Authorities

Without adjustment for deprivation After adjustment for deprivation Upper quartile Upper quartile 0.0

25.0

50.0

Rank 75.0

0.0

25.0

50.0

Rank 75.0

Posterior probability that HA is ranked in top 5%

Without adjustment (samples)means for p.bottom5

for deprivation and clinical stage N N After adjustment for deprivation (64 ) 0 .0 - 0 .1

and clinical stage (samples)means for p.bottom5

(6) 0.1 - 0.2

(4) 0.2 - 0.5

N 0.1-0.2

(14 ) < 0.0

(72 ) 0 .0 - 0 .1

(6) 0.1 - 0.2

(2) 0.2 - 0.5

(1) >= 0.5

200.0km

Outline

Joint spatial variation in risk of multiple diseases

Disease 1 Disease 2

SMR SMR 0.7

1.0

1.5

Specific component 1 SMR 0.7

1.0

1.5

Shared component Specific component 2 RR 0.7

1.0

1.5

Knorr-Held and Best (2001) 0.7

1.0

1.5

RR 0.7

1.0

1.5

Statistical model

Y 1i ~ Poisson( R 1i E 1i ); log R 1i = S i + U 1i Y 2i ~ Poisson( R 2i E 2i ); log R 2i = S i + U 2i S i ~ spatial model (shared component of risk) U 1i ~ spatial model (component of risk specific to disease 1) U 2i ~ spatial model (component of risk specific to disease 2) • Extends to >2 diseases (Tzala and Best, 2006) • Extends to shared variations in space and time (Richardson et al, 2006) © Imperial College London

Joint variation in COPD and lung cancer in GB

N values for SMR1 COPD SMR (42 ) < 0.66

(10 0) 0.66 - 0 .8

(13 9) 0.8 - 1 .0

(12 1) 1.0 - 1.25

(40 ) 1 .25 - 1.5

(17 ) >= 1.5

N values for SMR2 Lung cancer SMR (11 ) < 0.66

(10 1) 0.66 - 0 .8

(20 1) 0.8 - 1 .0

(10 3) 1.0 - 1.25

(33 ) 1 .25 - 1.5

(10 ) >= 1.5

200.0km

Best and Hansell (2009) 200.0km

Modelled risk estimates

values for shared N Shared risk (16 8) < 0.9

(40 ) 0 .9 - 0 .9 5 N values for log.COPD.spec.MALE.model2

COPD specific risk (80 ) < -0.1

(39 ) -0 .1 - -0.05

(71 ) -0.05 - 0 .0

(32 ) 0 .95 - 1.0

(11 1) 0.0 - 0 .0 5 (34 ) 1 .0 - 1 .0 5 Shared risk interpreted as (25 ) 1 .05 - 1.1

mainly reflecting (28 ) 1 .1 - 1 .1 5 geographical variations in (13 2) >= 1.15

community-level smoking behaviour COPD specific risk interpreted as reflecting smoking-adjusted variations in COPD mortality (89 ) 0.05 - 0 .1

(36 ) 0 .1 - 0.15

(33 ) >= 0.15

200.0km

MLE estimate of relative survival: breast N

Joint variation in relative survival of colon and breast cancer by English Health Authority

• Shared spatial patterns of relative survival may reflect variations in effectiveness of health care system N N Observed 5-year relative MLE estimate of relative survival: breast survival: Breast < 65% Observed 5-year relative (13) 65.0 - 70.0

survival: Colon N (36) 70.0 - 75.0

(38) 75.0 - 80.0

< 20% (5) >= 80.0

20% to 30% (2) < 20.0

(10) 20.0 - 30.0

(31) 30.0 - 40.0

(36) 40.0 - 50.0

(16) >= 50.0

200.0km

shared differential in relative survival N N

Difference in relative survival in each HA compared to England as a whole

Shared difference shared differential in relative survival N < -30% Posterior Prob that shared difference > 0 N (13) -3.0 - -1.5

< 0.2

0.2 – 0.8

(1) >= 3.0

> 0.8

(11) < 0.2

(73) 0.2 - 0.8

(11) >= 0.8

200.0km

breast specific differential in rel. surv N

Difference in relative survival in each HA compared to England as a whole

Difference specific to breast cancer breast specific differential in rel. surv N < -30% -15% to -30% N Difference specific to colon cancer (1) -3.0 - -1.5

N < -30% -15% to -30% -15% to 15% (9) < -3.0

(19) -3.0 - -1.5

(41) -1.5 - 1.5

(16) 1.5 - 3.0

(10) >= 3.0

200.0km

Diet-related cancers in Greece

N Spatial common factor values for p 1 N Cancer-specific spatial residuals oesophagus N (10) -0.1184 - -0.02204

stomach (8) -0.02204 - 0.03013

(8) 0.03013 - 0.1332

(12) >= 0.1332

(18) <-0.1163

(6) -0.1163 - -0.01743

(1) -0.01743 - 0.03017

(7) 0.03017 - 0.1402

(19) >= 0.1402

values for p N colorectal 200.0km

200.0km

N (10) -0.1234 - -0.01909

pancreas (7) -0.01909 - 0.0284

(17) 0.0284 - 0.1359

(3) >= 0.1359

(6) <-0.1234

(4) -0.1234 - -0.01909

(2) -0.01909 - 0.0284

(15) 0.0284 - 0.1359

(24) >= 0.1359

200.0km

Cut-points based on quintiles of distribution of factor values and of residuals across all cancers Tzala and Best (2006) values for p N prostate 200.0km

(30) -0.01909 - 0.0284

(3) 0.0284 - 0.1359

(0) >= 0.1359

bladder 200.0km

(10) <-0.1184

(9) -0.1184 - -0.02204

(14) -0.02204 - 0.03013

(15) 0.03013 - 0.1332

(3) >= 0.1332

200.0km

Outline

Extensions of disease mapping to space time modelling Noisy data in each area Noise model: Poisson/Binomial Latent structure: Space + Time + (Residuals)

+ joint Bayesian estimation

Inference

Basic space-time model set-up

Y it ~ Poisson(R it E it ); log R it = S i + T t + U it S i ~ spatial CAR model (common spatial pattern) T t ~ random walk (RW) model (common temporal trend) U it ~ Normal(0, v) (space-time residual reflecting idiosyncratic variation) • Extends to shared variations of 2 outcomes in space and time © Imperial College London

Space-time variations in Male and Female lung cancer incidence

(Richardson et al, 2006) • Lung cancer, with its low survival rates, is the biggest cancer killer in the UK – Over one fifth of all cancer deaths in UK are from lung cancer (25% for male and 18% for female) • Major risk factor is smoking. – Smoking time trends different for men/women: uptake of smoking started to decrease in cohorts of men after 1970, while for women the levelling off was later, after 1980 • Other risk factors include exposure to workplace agents, radon, air pollution …  Interested in similarity and specificity of patterns between men and women

Space-time analysis of Male and Female lung cancer incidence

Male/Female lung cancer incidence in Yorkshire:81-85, 86-90, 91-95, 96-99 (Richardson, Abellan, Best, 2006)

Shared and specific patterns and time trends

Shared component Female/Male differential Time trend for male RRs in 10 wards Time trend for female RRs in 10 wards

Detection of space-time interaction patterns

Detection of space-time interaction patterns Noisy data in each area Noise model: Poisson/Binomial Latent structure: Space + Time + (Residuals)

+ joint Bayesian estimation

Inference

Detection of space-time interaction patterns Noisy data in each area Noise model: Poisson/Binomial Latent structure: Space + Time + Interactions

Any patterns?

+ + joint Bayesian estimation

Inference

Detection of space-time interaction patterns

• Study the persistence of patterns over time – Interpreted as associated with stable risk factors, environmental effects, socio-economic determinants • Highlight unusual patterns, via the inclusion of space time interaction terms, which are modelled by a mixture model • Unusual patterns in some areas may be linked to recording changes, emerging environmental hazards, impact of new policy or intervention program, …  a general tool for surveillance ?

Detection of space-time interaction patterns

Y it ~ Poisson( R it E it ); log R it = S i + T t + U it S i ~ spatial CAR model (common spatial pattern) T t ~ random walk (RW) model (common temporal trend) U it ~ Normal(0, v) (space-time interaction; idiosyncratic variation) © Imperial College London

Detection of space-time interaction patterns

Y it ~ Poisson( R it E it ); log R it = S i + T t + U it S i ~ spatial CAR model (common spatial pattern) T t ~ random walk (RW) model (common temporal trend) U it ~ q Normal(0, v 1 ) +(1-q) Normal(0, v 2 ); v 2 > v 1 (mixture model to characterise ‘stable’ and ‘unstable’ patterns over time) • Compute posterior probability, p it , that interaction parameter U it comes from the Normal(0, v 2 ) component • Classify area as ‘unstable’ if p it > 0.5 for at least one time, t (simulation study → 10% false positive rate; 20% false negative rate) © Imperial College London

Detecting unusual trends in congenital anomalies rates in England (Abellan et al 2008)

• Annual postcoded data on congenital anomalies (non chromosomal) recorded in England for the period 1983 – 1998 • Annual postcoded data on total number of live births, still births and terminations • 136,000 congenital anomalies  84.5 per 10 5 • Congenital anomalies are sparse: birth-years  Grid of 970 grid squares with variable size, to equalize the number of births and expected cases per area • Variations could be linked to socio-economic or environmental risk factors or heterogeneity in recording practices  Interest in characterising space time patterns © Imperial College London

Congenital anomalies in England, 1983-1998

Spatial main effect: evidence of spatial heterogeneity, linked to deprivation and maternal age Temporal main effect: downward trend around 1990 reflects implementation of “minor anomalies” exclusion policy

Congenital anomalies: Space-time interactions

Most areas are stable ( cluster 1 ) Some have a change around 90-91 where modifications in the classification of anomalies occurred ( clusters 2 and 3 ) Identified one very unusual time profile due to a change of local recording practice

Outline

1. Introduction 2. Mapping and spatial smoothing of health data 3. Classifying areas and summarising geographical variations in health outcomes 4. Modelling and mapping multiple outcomes 5. Modelling and mapping temporal trends Summary 6. Hierarchical related regression models for combining individual and small area data

Summary

• Smoothing of small area risks is important to help separate ‘ signal ’ (spatial pattern) from ‘ noise ’ –Allows meaningful inference even when data are sparse • Achieved by ‘borrowing’ information neighbouring regions from • Bayesian hierarchical modelling provides formal method for carrying out this ‘borrowing of information’ –Provides rich output for statistical inference (estimation, quantification of uncertainty, hypothesis testing) –But, depends on “structural” assumptions built into the model (e.g. spatial dependence) –Computationally intensive © Imperial College London

Summary

Bayesian approach extends naturally to allow: • Adjustment for covariates (see later) • Joint mapping of 2 or more health outcomes • Joint modelling of spatial and temporal variation

Benefits of Joint Analysis of related health outcomes

Joint analysis of two related health outcomes is of interest in several contexts: • Epidemiology: quantify ‘expected’ variability linked to shared risk factors and tease out specific patterns • Health planning: assess the performance of the health system, e.g. for health outcomes linked to screening policies • Data quality issues: uncover anomalous patterns linked to a data source shared by several outcomes

Benefits of Space Time Analysis for (non-infectious) health outcomes

• Study the persistence of patterns over time – Interpreted as associated with stable risk factors, environmental effects, distribution of health care access … • Highlight unusual patterns in time profiles inclusion of space-time interaction terms via the – Time localised excesses linked to e.g. emerging environmental hazards with short latency – Variability in recording practices  Increased epidemiological interpretability  Potential tool for surveillance

Outline

Introduction

• Models and applications discussed so far have focused on: – describing outcomes geographical and temporal patterns in health – partitioning sources of variation into e.g. systematic and idiosyncratic, spatial and temporal, shared and specific, … • Growing interest in trying to explain variations at level of areas and geographical individuals → Build regression models and explanatory variables linking health outcomes

Regression models for small area data

Standard regression for individual-level outcomes Individual exposure

x ij y ij

Individual outcome Aggregate exposure

Z i , X i

Ecological regression Aggregate exposure

Z i , X i Y i

Aggregate outcome Hierarchical Related Regression (HRR; Jackson et al, 2006, 2008a,b) Individual exposure

x ij y ij

Individual outcome Aggregate exposure

Z i , X i Y i

Aggregate outcome

Case study: Socioeconomic inequalities in health

Jackson, Best and Richardson (2008b) • Geographical inequalities in health are well documented • One explanation is that people with similar characteristics cluster together, so area effects are just the result of differences in characteristics of people living in them ( compositional effect ) • But, evidence suggests that attributes of places may influence health over and above effects of individual risk factors ( contextual effect ) – economic, environmental, infrastructure, social capital/cohension Question • Is there evidence of contextual effects risk of limiting long term illness (LLTI) of area of residence on and heart disease, after adjusting for individual-level socio-demographic characteristics

Data and Methodological Issues

Methodological issues • Surveys typically contain sparse individual data difficult to estimate contextual effects per area so • Can’t separate individual and contextual effects using only aggregate data ( ecological bias ) • Improve power and reduce bias by combining data using new class of multilevel models developed by BIAS Our goal: data synthesis using • Individual-level data – Health Survey for England, 1997-2001 • Area-level (electoral ward) data – 1991 census small-area statistics – Hospital Episode Statistics

Data sources

INDIVIDUAL DATA Health Survey for England • Self-reported limiting long term illness • Self reported hospitalisation for heart disease • age and sex • ethnicity • social class • car access • income • etc.

Ward codes made available under special license AREA (WARD) DATA Census small area statistics • Carstairs deprivation index (area-level material deprivation) Individual-level Health outcomes Contextual effect Individual predictors

Multilevel model for individual data

b x ij y ij person j c

y ij = disease (1) / no disease (0) x ij = non-white (1) / white (0) Z i = deprivation score

m,v 2 a Z i i area i

y ij ~ Bernoulli(p ij ), person j, area i logit p ij a i = a i + b x ij + c Z i ~ Normal(m, v 2 ) b = relative risk of disease for non-white versus white individual c = contextual effects a i = “unexplained” area effects

Results from analysis of individual survey data: Heart Disease (n=5226)

0.2

0.5

1.0

2.0

odds ratio

Univariate regression

Area deprivation

Multiple regression

Area deprivation No car Social class IV/V Non white 5.0

Results from analysis of individual survey data: Limiting Long Term Illness (n=1155)

0.5

1.0

odds ratio 2

Univariate regression

Area deprivation

Multiple regression

Area deprivation Female Non white Doubled income 3

Comments

• CI wide and not significant for most effects • Some evidence of contextual effect of area deprivation for both heart disease and LLTI – Adjusting for individual risk factors (compositional effects) appears to explain contextual effect for heart disease – Unclear whether contextual effect remains for LLTI after adjustment for individual factors – Survey data lack power to provide reliable answers about contextual effects • What can we learn from aggregate data?

Contextual effect Aggregate health

Area-level data

AREA (WARD) DATA Census small area statistics • Carstairs deprivation index • population count by age and sex • proportion reporting LLTI • proportion non-white • proportion in social class IV/V • proportion with no car access PayCheck (CACI) • mean & variance of household income Hospital Episode Statistics • number of admissions for heart disease

Standard ecological regression model

Y i N i = number with disease = population X Z i i = proportion non-white = area deprivation score

M,V 2 B C

Y i ~ Binomial(q i , N i ), area i logit q i A i = A i + BX i + CZ i ~ Normal(M, V 2 ) B = association between disease prevalence and proportion non-white C = contextual effects A i = “unexplained” area effects

A i Z i area i Y i X i N i

Comparison of individual and ecological regressions: Heart Disease

Individual Ecological 0.2

0.5

1.0

2.0

odds ratio Area deprivation No car Social class IV/V Non white 5.0

Comparison of individual and ecological regressions: Limiting Long Term Illness

Individual Ecological Area deprivation Female Non white Doubled income 0.5

1.0

2 odds ratio 3 6

Ecological bias

Ecological bias (difference between individual and aggregate level effects) can be caused by: • Confounding – confounders can be area-level (between-area) or individual level (within-area).

→ include control variables and/or random effects in model • Non-linear covariate-outcome relationship, combined with within-area variability of covariate – No bias if covariate is constant in area (contextual effect) – Bias increases as within-area variability increases – …unless models are refined to account for this hidden variability

Standard ecological regression model

Y i N i = number with disease = population X Z i i = proportion non-white = area deprivation score

M,V 2 B C

A i Z i area i Y i X i N i

Integrated ecological regression model

Y i N i = number with disease = population Average of the X i Z i = area deprivation score ij , in area i

m,v 2

Y i q i a i ~ Binomial(q i , N i ), area i =  p ij (x ij ,Z i ,a i , b, c)f i (x)dx ~ Normal(m, v 2 ) b = relative risk of disease for non-white versus white individual c = contextual effects a i = “unexplained” area effects

a i Z i area i Y i b c X i N i

Combining individual and aggregate data

Multilevel model for individual data

b c m,v 2 x ij y ij person j a i Z i area i m,v 2

Integrated ecological model

b c a i Z i area i Y i X i N i

Combining individual and aggregate data

Hierarchical Related Regression (HRR) model (Jackson, Best, Richardson, 2006, 2008a,b)

b c m,v 2

Joint likelihood for y ij and Y i depending on shared parameters a i , b, c, m, v 2

x ij y ij person j a i Z i area i Y i X i N i

Combining individual and aggregate data

Hierarchical Related

b c m,v 2

Regression (HRR) model Joint likelihood for y ij and Y i depending on 2006, 2008a,b) shared parameters R software (maximum likelihood) 2 or WinBUGS (Bayesian)

a i x ij y ij person j Z i Y i X i N i area i

Comparison of results from different regression models: Heart Disease

Individual Standard ecological Integrated ecological HRR Area deprivation No car Social class IV/V PR95 = 10.1; 95% CI(5.3, 18.1) PR95 = 4.2; 95% CI(3.6, 5.1) 0.2

0.5

1.0

2.0

odds ratio 5.0

10 Non white

Comparison of results from different regression models: Limiting Long Term Illness

Individual Standard ecological Integrated ecological HRR PR75 = 2.7; 95% CI(1.7, 4.1) PR75 = 2.9; 95% CI(2.4, 3.7) Area deprivation Female 0.5

1.0

odds ratio 2 3 6 Non white Doubled income

Comments

• Integrated ecological model yields odds ratios that are consistent with individual level estimates from survey • Large gains in precision achieved by using aggregate data • Significant contextual effect of area deprivation for LLTI but not heart disease • More unexplained between-area variation (PR95) for heart disease than LLTI • Little difference between estimates based on aggregate data alone and combined individual + aggregate data – Individual sample size very small (~0.1% of population represented by aggregate data) • In other applications with larger individual sample sizes and/or less informative aggregate data , combined HRR model yields greater improvements (simulation study)

Strengths of HRR approach……

 Aims to provide individual-level inference using aggregate data by:  Fitting integrated individual-level model to alleviate one source of ecological bias  Including samples of individual data to help identify effects  Uses data from all geographic areas (wards, constituencies) , not just those in the survey  Improves precision of parameter estimates  Improves ability to investigate contextual effects

…..and limitations of HRR approach

 Integrated individual-level model relies on large contrasts in the predictor proportions across areas  e.g. limited variation in % non-white across constituencies: (median 2.7%, 95 th percentile 33)  Our estimates may not be completely free from ecological bias (Jackson et al, 2006)  If individual level data too sparse , may be overwhelmed by aggregate data

Data requirements for HRR models

• Individual data – requires geographical (group) identifiers for individual data • Aggregate data – requires large exposure contrasts between areas – requires information on within-area distribution covariates of • Important to check compatibility of different data sources when combining data

Thank you for your attention

Acknowledgements: Sylvia Richardson, Juanjo Abellan, Virgilio Gomez-Rubio, Chris Jackson Training courses in Bayesian Analysis of Small Area Data using WinBUGS and INLA, London, July 13-16 2010 See www.bias-project.org.uk

for details

References Abellan JJ, Richardson S and Best N. Use of space-time models to investigate the stability of patterns of disease. Environ Health Perspect 116(8), (2008), 1111-1119.

Best N and Hansel A. Geographic variations in risk: adjusting for unmeasured confounders through joint modelling of multiple diseases. Epidemiology, 20(3), (2009), 400-410. Best, N.G., Richardson, S. and Thomson, A. A comparison of Bayesian spatial models for disease mapping. Statistical Methods in Medical Research 14, (2005), 35-59.

Jackson C, Best N and Richardson S. Hierarchical related regression for combining aggregate and individual data in studies of socio-economic disease risk factors. J Royal Statistical Society Series A: Statistics in Society 171 (2008a) ,159-178 Jackson C, Best N and Richardson S. Studying place effects on health by synthesising individual and area level outcomes. Social Science and Medicine, 67, (2008b), 1995-2006 Jackson C, Best N and Richardson S. Improving ecological inference using individual-level data. Statistics in Medicine, 25, (2006), 2136-2159 Knorr-Held, L. and Best, N.G., A shared component model for detecting joint and selective clustering of two diseases, Journal of the Royal Statistical Society, Series A 164, (2001), 73-86.

Richardson, S., Thomson, A., Best, N.G. and Elliott, P. Interpreting posterior relative risk estimates in disease mapping studies. Environmental Health Perspectives 112, (2004), 1016-1025.

Richardson, S., Abelan, J.J.,and Best, N. Bayseian spatio-temporal analysis of joint patterns of male and female lung cancer in Yorkshire (UK). Statistical Methods in Medical Research 15, (2006), 385-407.

Tzala, E. Best, N. Bayesian latent variable modelling of multivariate spatio-temporal variation in cancer mortality. Stat Methods Med Res 17, (2008), 97-118.

Spanish disease atlases and other related resources López-Abente G, Pollan M, Escolar A, Errezola M, Abraira V. 1996. Atlas of cancer mortality and other causes of death in Spain 1978 1992. Madrid: Fundación Científica de la Asociación Española contra el Cáncer. http://www2.uca.es/hospital/atlas92/www/Atlas92.html

. Maps of SMRs and Directly Standardised Rates at province level Benach J, Yasui Y, Borrell C, Rosa E, Pasarin MI, Benach N, Español M, Martinez JM, Daponte A. 2001. Atlas of mortality at small area level in Spain 1987-1995. Barcelona: Universitat Pompeu Fabra. Empirical Bayes smoothing of relative risks at municipality level (aggregated some municipalities to reduced sparseness).

All of the following produce maps of Bayesian smoothed relative risks at municipality level: López-Abente G, Ramis R, Pollán M, Aragonés N, Pérez-Gómez B, Gómez-Barroso D, Carrasco JM, Lope V, García-Pérez J, Boldo E, García-Mendizabal MJ. 2007. Atlas of cancer mortality at municipality level in Spain 1989-1998. Área de Epidemiología Ambiental y Cáncer del Centro Nacional de Epidemiología, ISCIII. Botella P, Zurriaga O, Posada de la Paz M, Martinez-Beneito MA, Bel E, Robustillo A, Ramalle E, Duran E, Sanchez-Porro P. 2006. National-Provincial Atlas of Rare Diseases 1999-2003. Martinez-Beneito MA, Lopez-Quilez A, Amador A, Melchor I, Botella P, Abellan C, Abellan JJ, Verdejo F, Zurriaga O, Vanaclocha H, Escolano M. 2005. Atlas of Mortality in the Valencian Region 1991-2000. Benach J, Matinez JM, Yasui Y, Borrell C, Pasarin MI, Español E, Benach N. 2004. Atlas of mortality at small area level in Catalonia 1984-1998. Barcelona: Universitat Pompeu Fabra / Fundació Jaume Bofill / Editorial Mediterrània

Spanish disease atlases and other related resources DEMAP group, Andalusian School of Public Health, Granada. Produced interactive mortality atlas for Andalucia, and socioeconomic indices at municipality level for Spain. See http://www.demap.es/Demap/index.html

MEDEA: Research network on Epidemiology and Public Health, working on socioeconomic and environmental inequalities in health at small area level. See http://www.proyectomedea.org/medea.html

VPM Atlas Project ( Atlas de Variaciones en la Práctica Médica). Studying and mapping variations in provision and usage of health care at small area level in 16 of the 17 regions of Spain, using data on hospital discharges. See http://www.atlasvpm.org/avpm/inicio.inicio.do

Smoothing of the RRs of hot spots (4 contiguous areas with average expected counts ≈ 5) for different spatial models Richardson et al (EHP, 2004) True RR = 3 True RR = 2