Introduction to Spatial Regression

1. Introduction to Spatial Regression Glen Johnson, PhD Lehman College / CUNY School of Public Health glen.johnson@lehman.cuny.edu

2. Typical scenario: • Have a health outcome and covariables aggregated at a common geographic level, such as counties, census tracts, ZIP codes … • Want to measure association between the outcome and the covariables. • Specific Question is: Are there variables that co-vary spatially with the outcome variable ?

3. Benzene in ambient air Smoking rate + Lung Cancer Rates (observed) + … + ? = + residual

4. Consider the linear model: This is the point of departure.

5. When applying regression modeling to spatial units that are connected in space (lattice data), the critical assumption that residuals are independently distributed with constant variance is typically violated. Tobler’s First Law of Geography: Things closer in space tend to be more similar than things further apart

6. When we model the expected value, E[y], as a function of spatially-varying covariates, it is possible that we may explain all of the spatial variation of the observed response, y, with the covariates, leaving uncorrelated residuals. When this is not the case, as is typical, the assumption of iid residuals is violated and we will obtained biased estimates of the variance – typically biased downward, leading to underestimating our standard errors and concluding that some covariates are significant when in fact they are not.

7. Tests for spatial autocorrelation should be applied to residuals if a “conventional” regression model is applied. This may be done with various software packages or GIS add-ons. A common statistic is Moran’s I, which equals

8. When residual spatial autocorrelation is present, several approaches may be taken to adjust for it. The simplest is to add a fixed effect dummy variable to allow the model intercept to change with spatial location. For example, an adjustment is made that depends on a county of membership. This is essentially stratifying the analysis by location And can be done with any statistical software

9. Since spatial location is a proxy for unobserved randomly varying covariables, it is more correctly treated as a random effect in a mixed effect model, such as Which can be solved for through pseudo-likelihood methods, using software likePROC GLIMMIX or PROC MIXED in SAS, or R with appropriate library (?)

11. Illustration: Community Teen Pregnancy Rates vs. Socioeconomic Status and Demographic Composition

12. For each ZIP code: • Response (i.e. Teen Pregnancy cases) • Predictors: • % pop. > age 24 w/ 4-year or greater college degree • % single-parent households out of households w/ at least one child < 18 years old • % of tot. pop. that is Black Alone • % of tot. pop. that is Hispanic, regardless of race • % of tot. pop. that is a foreign-born naturalized citizen • % of tot. pop. with income below poverty • Population at Risk • County(crude indicator of neighborhood effect)

13. . . .

15. The Model … For i = 1, …,n ZIP codes, let yi= observed caseload ni= population at risk {x1, …, xp}i= community predictors {β1, …, βp} = coefficients Li = location effect, arising from a random process such that Li~ N(0, σL2) Then, the expected value of yi, given {x1, …, xp, L}i = E[yi| {x1, …, xp, L}i ] = niexp(β1x1i+ … + βp xpi + Li)

16. Values for the unknown coefficients {β1, …, βp, σL2} are estimated with SAS PROC GLIMMIX, assuming yiarose from a Poisson random process, conditional on location. • … thus allowing risk adjusted estimates of caseload for each ZIP code. • Incorporating the “location effect”- adjusts for unidentified covariables that co-vary spatially with the response, thus reducing residual spatial autocorrelation and potential confounding- also provides a “smoothing” effect, in that the predicted caseload is adjusted towards a common local value

24. Other approaches include … A spatial lag model, where and a spatial error model, where for a spatial autoregressive coefficient ρ. These two models differ by whether the adjustment is made by a weighted sum of the response variable or the residuals.

25. The spatial lag and spatial error models can be solved for in Geoda, a simple, well supported freeware found at http://geodacenter.asu.edu/ … but only for gaussian responses. For generalized linear models (i.e. Poisson and logistic regression), see R with appropriate libraries http://www.r-project.org/

26. Another approach is hierarchical modelling, which treats the response as conditional on the weighted average of local neighborhood errors.

27. Frequentist solutions exist, but these hierarchical models lend themselves well to a fully Bayesian solution, as used by many geographic epidemiologists Main advantages include * flexibility offered by Generalized Linear Mixed Models * obtain full distribution of possible outcomes - allows many ways to view the outcome (mean, median, percentiles) - inference based on actual probability distributions, instead of confidence intervals Main limitation is level of conceptual difficulty; however, implementation is accessible through free software … WINBUGS (Bayesian Inference Using the Gibbs Sampler)

28. A Hierarchical Model

29. Focus on the random effect that captures • local spatial autocorrelation • is distributed conditionally on location, such that

30. A Directed Acyclic Graph of the Bayesian Model

31. Gibbs sampling basic procedure • All stochastic parameters in the model are assigned an initial value (somewhat arbitrarily). • The values for each parameter are updated by random simulation from a conditional probability distribution, given all other parameters in the model. • After all terms have been updated, completing one cycle (of what is called a Markov Chain), the cycle is repeated. • After many iterations, the simulated values for each term converge to a stationary posterior distribution (further iterations don’t change the distribution) • Estimation and inference can then be made from these posterior distributions • For example, a simulated sample of 1000 fitted SIR values (μi / Ei) can be used to yield a point estimate (typically the median)and an interval estimate, such as the 95 %-tile range (credible set)

32. 5th %-tile 50th %-tile 95th %-tile

33. An illustration for geospatial analysis of prostate cancer incidence in New York State, USA …

34. Prostate Cancer Incidence by ZIP codeadjusted for age and race New York State1994-1998

35. Example Output: Posterior Kernel Densities of Prostate Cancer Incidence (`94-`98) for Some Manhattan ZIP Codes

39. some references • Waller, L.A. and Gotway, C.A. 2004. Applied Spatial Statistics for Public Health Data. Wiley. 494 pp. • Johnson, G.D. 2004. Smoothing Small Area Maps of Prostate Cancer Incidence in New York State (USA) using Fully Bayesian Hierarchical Modelling. International Journal of Health Geographics 2004, 3:29 ( http://www.ij-healthgeographics.com/content/3/1/29 ) • Elliot, P., Wakefield, J.C., Best, N.G. and Briggs, D.J. 2000. Spatial Epidemiology: Methods and Applications. Oxford. 475 pp. • Statistics in Medicine. 2000. Vol. 19 (special issue on disease mapping) • Lawson, A. et al. 1999. Disease Mapping and Risk Assessment for Public Health. Wiley. 482 pp.

40. Method and Software Sources GeoDa http://geodacenter.asu.edu/ (with links to R and R-Geo) WINBUGS for Bayesian Modeling http://www.mrc-bsu.cam.ac.uk/bugs/welcome.shtml Both of these freewares are supported by large international community with active listserves