Spatial Modeling of Trends in Crime over Time in Philadelphia

TL;DR

This paper develops Bayesian spatial models to analyze and predict crime trends over time in Philadelphia, incorporating neighborhood predictors and spatial discontinuities for improved accuracy.

Contribution

It introduces novel Bayesian approaches for spatial sharing of information in crime trend modeling, including handling spatial discontinuities.

Findings

Conditional autoregressive models improve predictive accuracy.

Spatial discontinuities may indicate natural barriers or built environment effects.

Incorporating economic and demographic predictors enhances model insights.

Abstract

Understanding the relationship between change in crime over time and the geography of urban areas is an important problem for urban planning. Accurate estimation of changing crime rates throughout a city would aid law enforcement as well as enable studies of the association between crime and the built environment. Bayesian modeling is a promising direction since areal data require principled sharing of information to address spatial autocorrelation between proximal neighborhoods. We develop several Bayesian approaches to spatial sharing of information between neighborhoods while modeling trends in crime counts over time. We apply our methodology to estimate changes in crime throughout Philadelphia over the 2006-15 period, while also incorporating spatially-varying economic and demographic predictors. We find that the local shrinkage imposed by a conditional autoregressive model has substantial benefits in terms of out-of-sample predictive accuracy of crime. We also explore the possibility of spatial discontinuities between neighborhoods that could represent natural barriers or aspects of the built environment.

Figures

Tables

Table 1. Table 1 : Comparison of predictive accuracy between the different models outlined in Section  3 . The mean squared error for both in-sample and out-of-sample predictions are provided, as well as the percentage change in MSE out subscript MSE out {\rm MSE}_{\rm out} relative to model ( 3.5 ) without shrinkage. We also provide the Moran’s I 𝐼 I measure of spatial correlation calculated on the estimated time trends β i subscript 𝛽 𝑖 \beta_{i} from each model.

Model	${\mathrm{MSE}}_{\mathrm{in}}$	${\mathrm{MSE}}_{\mathrm{out}}$	$\%$ change	${\mathrm{MSE}}_{\mathrm{cv}}$	Moran’s $I$
Global $\alpha ,\beta$	0.3558	0.3694	+182.4	0.3043	-
Separate ${\alpha }_i,{\beta }_i$ Models
No Shrinkage (3.5)	0.0567	0.1308	-	0.1001	0.17
Global Shrinkage (3.6)-(3.7)	0.0698	0.1080	-17.4	0.0928	0.17
Spatial CAR (3.11)-(3.12)	0.0703	0.1052	- 19.5	0.0922	0.61
Variable Borders (3.13)-(3.14)	0.0706	0.1069	-18.2	0.0927	0.49

Table 2. Table S1 : Comparison of predictive accuracy between the different models outlined in Section 3 with non-informative priors on the hyper parameters. The mean squared error for both in-sample and out-of-sample predictions are provided. We also provide the Moran’s I 𝐼 I measure of spatial correlation calculated on the estimated time trends β i subscript 𝛽 𝑖 \beta_{i} from each model.

Model	${\mathrm{MSE}}_{\mathrm{in}}$	${\mathrm{MSE}}_{\mathrm{out}}$	${\mathrm{MSE}}_{\mathrm{cv}}$	Moran’s $I$
Separate ${\alpha }_i,{\beta }_i$ Models
Global Shrinkage	0.0698	0.1080	0.0927	0.17
Spatial CAR	0.0701	0.1052	0.0922	0.61
Variable Borders	0.0706	0.1069	0.0927	0.48

Table 3. Table S2 : Estimate and standard error for each coefficient γ j subscript 𝛾 𝑗 \gamma_{j} . For the Bayesian models, we report the mean and the standard deviation from 1000 independent draws from the posterior distribution.

	No shrinkage		Global shrinkage		Spatial CAR		Variable borders
	Estimate	St.Error	Mean	SD	Mean	SD	Mean	SD
log.income	-0.186	0.009	-0.186	0.022	-0.113	0.023	-0.099	0.021
sqrt.poverty	0.182	0.009	0.182	0.023	0.095	0.020	0.110	0.018
segregation	0.010	0.005	0.010	0.013	-0.023	0.017	-0.013	0.017
sqrt.vacantprop	0.116	0.006	0.115	0.014	0.052	0.015	0.045	0.014
sqrt.comresprop	0.227	0.005	0.223	0.013	0.241	0.011	0.240	0.011
pop.total	0.216	0.005	0.212	0.013	0.263	0.011	0.317	0.012

Equations

$${\text{segregation}}_i=\frac{1}{2}\sum _r|p_{i,r}-{\overline{p}}_r|$$

$${\text{poverty}}_i=\sum _{j=1}^7{w}_j{q}_{i,j}$$

$${\text{vacancy}}_i=\frac{Are{a}_i(\text{Vacant})}{Are{a}_i}$$

$${\text{comresprop}}_i=\frac{Are{a}_i(\text{Commercial})}{Are{a}_i(\text{Commercial})+Are{a}_i(\text{Residential})}$$

$$I=\frac{n}{{\sum }_i{\sum }_j{w}_{ij}}\frac{{\sum }_i{\sum }_j{w}_{ij}(X_i-\bar{X})(X_j-\bar{X})}{{\sum }_i(X_i-\bar{X}{)}^2}$$

$$y_{it}=\log (c_{it}+\sqrt{c_{it}^2+1})-\log (2)$$

$$y_{it}=\alpha +{\mathbf{z}}_i^T\mathbf{\gamma }+e_{it},$$

$$y_{it}=\alpha +{\mathbf{z}}_i^T\mathbf{\gamma }+\beta \cdot t+{ϵ}_{it}\mathrm{where}{ϵ}_{it}\sim \mathrm{N}(0,{\sigma }^2)$$

$$y_{it}={\alpha }_i+{\mathbf{z}}_i^T\gamma +{\beta }_i\cdot t+{ϵ}_{it}\mathrm{where}{ϵ}_{it}\sim \mathrm{N}(0,{\sigma }^2).$$

$$y_{it}={\alpha }_i+{\beta }_i\cdot t+{ϵ}_{it}\mathrm{where}{ϵ}_{it}\sim \mathrm{N}(0,{\sigma }^2)$$

$$\mathbf{\alpha }\sim \mathrm{N}({\alpha }_0\cdot 1,{\tau }_{\alpha }^2\cdot \mathbf{I})$$

$$\mathbf{\beta }\sim \mathrm{N}({\beta }_0\cdot 1,{\tau }_{\beta }^2\cdot \mathbf{I})$$

$$\mathbf{\gamma }\sim \mathrm{N}(0,{\tau }_{\gamma }^2\cdot \mathbf{I})$$

$${\sigma }^2$$

$${\tau }_{\alpha }^2$$

$${\tau }_{\beta }^2$$

$${\tau }_{\gamma }^2$$

$${\theta }_i\mid {\mathbf{\theta }}_{-i},{\theta }_0,{\tau }^2\sim \mathrm{N}\left(\frac{\rho {\sum }_j{w}_{ij}{\theta }_j+(1-\rho ){\theta }_0}{\rho {\sum }_j{w}_{ij}+(1-\rho )},\frac{{\tau }^2}{\rho {\sum }_j{w}_{ij}+(1-\rho )}\right),$$

$$\mathbf{\theta }|{\theta }_0,{\tau }^2\sim \mathrm{N}\left({\theta }_0\cdot 1,{\tau }^2\cdot [\rho ({\mathbf{D}}_W-\mathbf{W})+(1-\rho )\mathbf{I}{]}^{-1}\right)$$

$$\mathbf{\alpha }\sim \mathrm{N}\left({\alpha }_0\cdot 1,{\tau }_{\alpha }^2\cdot \mathbf{\Sigma }\right)$$

$$\mathbf{\beta }\sim \mathrm{N}\left({\beta }_0\cdot 1,{\tau }_{\beta }^2\cdot \mathbf{\Sigma }\right)$$

$$\mathbf{\alpha }|{\mathbf{W}}^{\alpha }\sim \mathrm{N}\left({\alpha }_0\cdot 1,{\tau }_{\alpha }^2\cdot {\mathbf{\Sigma }}_{\alpha }\right)$$

$$\mathbf{\beta }|{\mathbf{W}}^{\beta }\sim \mathrm{N}\left({\beta }_0\cdot 1,{\tau }_{\beta }^2\cdot {\mathbf{\Sigma }}_{\beta }\right)$$

$${\mathrm{MSE}}_{\mathrm{in}}=\frac{1}{1336}\sum _{i=1}^{1336}(y_{i,2014}-{\hat{y}}_{i,2014}{)}^2.$$

$${\mathrm{MSE}}_{\mathrm{out}}=\frac{1}{1336}\sum _{i=1}^{1336}(y_{i,2015}-{\hat{y}}_{i,2015}{)}^2.$$

$${\mathrm{MSE}}_{\mathrm{cv}}=\frac{1}{10}\sum _{t=1}^{10}{\mathrm{MSE}}_{\mathrm{out}}^t\mathrm{where}{\mathrm{MSE}}_{\mathrm{out}}^t=\frac{1}{1336}\sum _{i=1}^{1336}(y_{i,t}-{\hat{y}}_{i,t}{)}^2.$$

$$\mathbf{\theta }|\mathbf{Y},{\mathbf{\theta }}_0,{\sigma }^2,{\mathbf{\tau }}^2\sim N(\hat{\mathbf{\theta }},{\mathbf{V}}_{\mathbf{\theta }})$$

$$\hat{\mathbf{\theta }}$$

$${\mathbf{V}}_{\mathbf{\theta }}$$

$${\alpha }_0|\mathbf{\alpha },{\tau }_{\alpha }^2$$

$${\beta }_0|\mathbf{\beta },{\tau }_{\beta }^2$$

$${\sigma }^2$$

$${\tau }_{\alpha }^2$$

$${\tau }_{\beta }^2$$

$${\tau }_{\gamma }^2$$

$${\sigma }^2|\mathbf{Y},\mathbf{\theta }$$

Full text

urbanism; crime; spatial; time trends,

1901.08117 \startlocaldefs

\endlocaldefs

and

1 Introduction

Modeling and prediction of crime has always been of interest to local authorities, police departments and governments to assure safety of the population and more efficient law enforcement. Recent availability of detailed crime data has made this effort even more accessible to statistical practitioners and the general public.

As an example, the Philadelphia police department has released detailed information about reported crimes committed from 2006 to the present day111http://www.phlcrimemapper.com/. The information about each reported crime includes the type of crime (which we will describe in Section 2), the date and time of the crime and the GPS location of the crime.

Using their reported crime data, many police departments have used statistical modeling procedures and algorithms to help predict locations of crimes for better prevention and faster intervention (Hvistendahl, 2016). The modeling of crime locations is not only useful for law enforcement but also for marketing strategies related both to real estate and commercial activities, e.g. Trulia222https://www.trulia.com/ uses crime data as part of their evaluation of the relative safety and attractiveness of different neighborhoods. In this paper, we will focus on estimating changes in violent crimes over the past decade at a local neighborhood resolution which will involve both temporal and spatial modeling of crime.

Many different approaches have been taken to the modeling of the spatial distribution of crime. These approaches can be subdivided into two general categories, either modeling crime as a spatial point process using the specific locations of each reported crime (Mohler et al., 2011; Taddy, 2010; Flaxman, 2014) or modeling crime as areal data, i.e. totals aggregated within larger regions, as in Aldor-Noiman et al. (2016), Law, Quick and Chan (2014) and Li et al. (2014).

A common method for modeling spatial point processes is kriging or Gaussian process interpolation (Stein, 2012; Cressie, 1990). This can be studied either with a classical approach, or with a Bayesian approach (Banerjee, Carlin and Gelfand, 2014). Alternative popular models consider other frameworks such as Gibbs point processes, Poisson processes and Cox processes; see Møller and Waagepetersen (2007).

Common classical methods for modeling areal data are spatial autoregressive models, that include the Simultaneous Autoregressive Model (Whittle, 1954), the spatial Durbin model (Anselin, 1998) and the Conditionally Auto Regressive model (Besag, 1974); for a review of these and other methods see LeSage and Pace (2009). Many of these models have also been considered and used in a Bayesian framework (Banerjee, Carlin and Gelfand, 2014).

Our goal in this paper is the estimation of trends in violent crime over the past decade at a high resolution local neighborhood level throughout the city of Philadelphia. As it is well established that crime frequencies are spatially correlated (Herbert, 1982; Brantingham and Brantingham, 1984), we need to create a model that allows the change in crime over time to be correlated by locally proximal neighborhoods. Our model will also account for characteristics of each local neighborhood, including the population count of the area and economic health of residents, as measured by median income and poverty level of households.

In addition to aiding law enforcement, accurate estimation of changes in crime at the local neighborhood level would also enable the study of the association between crime trends and changes in the built environment. We are particularly interested in how aspects of the built environment encourage vibrancy, a measure of positive human activity, and how vibrancy is associated with safety in local neighborhoods (Humphrey et al., 2017).

The city of Philadelphia is a particularly interesting case study for estimating trends in crime as it is a large urban area that is currently undergoing substantial development and experiencing population growth for the first time in decades. In addition to our primary goal of estimation of changes in crime in Philadelphia neighborhoods, this application also provides an interesting spatio-temporal data context for comparing different Bayesian shrinkage approaches to spatial areal modeling.

We will take an areal approach to modeling crime since our primary goal is greater understanding of evolving crime dynamics at the local neighborhood level within the city of Philadelphia. Our areal units will be U.S. Census block groups which consist of 10-20 city blocks and which are naturally interpretable as “neighborhoods”. U.S. Census block groups are also the highest resolution for which economic data is available as covariate information.

Compared to previous areal approaches (e.g. Aldor-Noiman et al. (2016), Law, Quick and Chan (2014) and Li et al. (2014)), we are using smaller areal units and we will focus on not only total crime but also the trend in crime over time within each local neighborhood. We have a longer time period (ten years) of recorded crimes for estimating time trends than Law, Quick and Chan (2014) that worked with property crimes over a two year period.

Our methodological contribution is the development of a Bayesian spatial modeling framework to explore global vs. local smoothing for our parameter estimates while also allowing for data-driven discontinuities in our model between proximal areal units. Using a Bayesian approach allows us to induce this smoothing through shrinkage priors for our parameters and also enables us to estimate borders between neighborhoods that have a high probability of being barriers.

In Section 2, we provide details for the neighborhood structure of Philadelphia and describe the detailed crime data that we will use to estimate changes in crime over the past decade. We also outline the demographic, economic and land use measures we will use as neighborhood-level predictors of violent crime in our spatial models. The code for acquiring and cleaning the data that were used in this analysis is available as a GitHub repository at https://github.com/cecilia-balocchi/Urban-project. In Section 3, we develop several Bayesian modeling approaches for global or local sharing of information between Philadelphia neighborhoods, as well as a model extension that allows for spatial discontinuities in our parameter estimates between proximal neighborhoods. We then compare these modeling options in terms of both in-sample and out-of-sample predictive accuracy in Section 4. We visualize and discuss the results of our spatial modeling of crime trends for Philadelphia in Section 5 and then conclude with a brief discussion in Section 6.

2 Population, Economic and Crime Data in Philadelphia

The population and economic data are provided by the US Census Bureau whereas crime data is provided by the Philadelphia Police Department. Our definition of local neighborhoods in Philadelphia will be based upon the “block group” geographical units defined by the US Census Bureau. The city of Philadelphia is divided into 384 census tracts which are divided into 1336 block groups. Shapefiles from the US Census Bureau give the boundaries and area of each census block group. Figure S1 in our supplementary materials (Balocchi and Jensen, 2019) gives a map outlining the 1336 block groups in Philadelphia.

Our motivation for analyzing trends in crime at this resolution is two-fold: a. US census block groups consist of 10-20 city blocks which generally matches our concept of a “neighborhood” and b. the block group level is the highest resolution of the economic data that we will use as predictors of crime. The average size of block groups in Philadelphia is 0.26 km2, with an average population of 1142 residents.

Our population data was pulled from the census website333https://factfinder.census.gov/ by setting the geography as all blocks in Philadelphia and setting the data source as “Hispanic or Latino Origin By Race” (which is SF1 P5 in their database). The raw demographic data gives the population count in each block group from the 2010 census. Figure S1 in our supplementary materials (Balocchi and Jensen, 2019) gives the population count for each block group in Philadelphia.

The same data also has the population count in each block group divided by ethnic categories444The ethnic categories are: White, Black, Asian, Native Americans, Native Pacific Islanders (including Hawaii), Other, Two or more races (nonhispanic) and Hispanic/Latino. We combined Native Americans, Native Pacific Islanders, and Two or more races into the Other category, which leads to five ethnicities in our analysis: 1. White, 2. Black, 3. Hispanic, 4. Asian, and 5. Other.. From these ethnicity counts, we calculate a measure of the segregation in each block group as

[TABLE]

where $p_{i,r}$ is the proportion of ethnicity $r$ in block group $i$ and $p_{r}$ is the proportion of ethnicity $r$ across the entire city of Philadelphia. The fraction $\tfrac{1}{2}$ scales this segregation measure to be between 0 and 1.

In addition to population count and our segregation measure, we will also consider several measures of the economic health of each neighborhood. Our economic data comes from the American Community Survey from the same US census website as our population data, specifically tables B19301 for income and C17002 for poverty, both from 2013. This data is only available at the resolution of census block groups. For each block group (neighborhood) in Philadelphia, we have income per capita as one predictor of crime.

We also have information about the proportion of households in various states of poverty. Specifically, we have the fraction of the population in seven different brackets of income-to-poverty-line ratios: $[0,0.5)$, $[0.5,1)$, $[1,1.25)$, $[1.25,1.5)$, $[1.5,1.85)$, $[1.85,2)$, $[2,\infty)$. For example, the $[0.5,1)$ bracket represents families that have income between $50\%$ of the poverty line and the poverty line itself. The poverty line is defined by the Census Bureau according to the size and composition of a household (e.g. a family with two children has a poverty line threshold of $23,999).

We use this poverty data to create a single measure of poverty for each block group (neighborhood) by calculating a weighted sum of the proportion of households in each of the seven poverty brackets:

[TABLE]

where $q_{i,1}$ is the proportion of households in block group $i$ that are in the lowest bracket $[0,0.5)$ and $q_{i,7}$ is the proportion of households in block group $i$ in the highest bracket $[2,\infty)$. We use linearly decreasing weights $\mathbf{w}=[1,5/6,4/6,3/6,2/6,1/6,0]$ to give higher weight to the brackets with higher poverty. Our poverty measure varies from 0 to 1, with larger values implying higher poverty.

In addition to the demographic and economic predictors described above, we also derive measures of the built environment that may also be predictive of crime. Our data on the built environment comes from the zoning designation of each lot in Philadelphia. Zoning data from the City of Philadelphia provides the area and registered land use designation (e.g. commercial, residential, industrial, vacant, transportation, park, civic) of all 560,000 lots in Philadelphia.

We create several land use metrics from these zoning designations that could be predictive of crime. First, we calculate the fraction of area in each block group $i$ that is designated as ‘Vacant’,

[TABLE]

Second, we calculate the ratio of the area in each block group $i$ that is commercial versus residential,

[TABLE]

To summarize, we have created six neighborhood characteristics that we will use as predictors of crime: population count, segregation, median household income, poverty, vacant proportion and commercial vs. residential proportion. Some block groups in Philadelphia have missing values for the economic predictors due to a very small or zero population count. We exclude these block groups (a total of eight) from our analysis. We additionally exclude one block group containing the detention centers in Philadelphia.

Our crime data comes from the Philadelphia Police Department and includes all crimes reported by the police in the city of Philadelphia from January 1, 2006 to December 31, 2015. For each reported crime, we have the type of crime, the date and time of the crime, and the location of the crime in terms of the GPS latitude and longitude (WGS84 decimal degrees). Each crime in our dataset is categorized into one of several types: homicide, sex crime, armed robbery, assault, burglary, theft, motor vehicle theft, etc.

We make a distinction between violent and non-violent (property) crimes in our analysis. As defined by the Uniform Crime Reporting program of the FBI , violent crimes include homicides, rapes, robberies and aggravated assaults whereas non-violent crimes include burglaries, thefts and motor vehicle thefts.

Our own crime categorization differs from the FBI in two ways. We combine ‘rapes’ and ‘sex assaults’ (which changed in definition in 2013) into a broader ‘sex crimes’ category and consider all ‘sex crimes’ as violent crimes. The FBI also makes a distinction between ‘aggravated assaults’ and ‘other assaults’, with the latter being where an injury does not occur but the threat of injury is present. In contrast, we combine both ‘aggravated assaults’ and ‘other assaults’ into a broader ‘assaults’ category and consider all ‘assaults’ as violent crimes.

For this paper, we focus entirely on the modeling of violent crimes as they have the most direct impact on human safety and the perception of safety. However, non-violent crimes are also important to track for law enforcement and are a focus of ongoing research. In the subsequent analyses in this paper, we will use ‘crime’ to mean only violent crimes.

In Figure 1, we give the counts of each type of violent crime within each year in 2006-2015, aggregated over the entire city. We see generally decreasing trends within the assault and robbery categories, which are the most numerous types of crimes. Sex crimes and homicides are also somewhat decreasing over this time span though it is harder to see this trend given the low counts for either type of crime.

Clearly, the impression given from Figure 1 is that violent crimes are generally decreasing in the city of Philadelphia over the time period from 2006 to 2015. However, are there specific neighborhoods that show substantially larger decreases or even some neighborhoods that show increases in violent crimes in this period?

As discussed in Section 1, we will model the spatial distribution of crime with an areal approach where our areal units are U.S. Census block groups which we define as the local neighborhoods of Philadelphia. Violent crimes are aggregated within each U.S. Census block group based on the GPS coordinates of each reported crime.

One issue with this approach is that some crimes occurring near to a boundary between U.S. Census block groups could be aggregated into the incorrect areal unit due to measurement error or ambiguity in their recorded point locations. This possibility is one of several motivations for our hierarchical Bayesian modeling approach that shares information between adjacent block groups when estimating crime totals and trends in crime over time across the city of Philadelphia.

In Figure S2 of our supplementary materials (Balocchi and Jensen, 2019), we give the count of violent crimes per year in each block group averaged over the years 2006-2015. One can see substantial heterogeneity across block groups in the average counts of violent crimes per year. There are several outlying values: particular block groups that have much higher average violent crime counts.

These outlying neighborhoods motivate us to examine violent crime totals on the log scale. In Figure S2 of our supplementary materials (Balocchi and Jensen, 2019), we also give the average of the logarithm of the count of violent crimes per year in each block group, averaged over the years 2006-2015. We can see more details of the spatial distribution of violent crime on the log scale. Modeling crime on the log scale has the additional benefit that changes in log crime can be interpreted as percentage changes in crime.

We also see in Figure S2 evidence of spatial correlation in violent crime totals between proximal block groups throughout the city. This is not surprising since the factors that lead to crime likely vary throughout the city in a (mostly) spatially continuous fashion. It is this spatial correlation that will be the focus of our modeling work in Section 3.

To get an idea of the strength of this spatial correlation, one of the standard statistics used for areal data is Moran’s I (Moran, 1950; Banerjee, Carlin and Gelfand, 2014), which is defined as

[TABLE]

where $W=(w_{ij})$ is a matrix of weights that capture the spatial proximity of the areal regions. We set $w_{ij}$ to be 1 if block groups $i$ and $j$ share a border and [math] otherwise. We use the queen contiguity method so two block groups share a border if they share at least a point on their boundaries.

Moran’s I can be used for testing for spatial autocorrelation: under the null hypothesis of no spatial association, we can compute exactly the mean (equal to $-\frac{1}{n-1}$) and standard error of Moran’s I. Calculating I on the total number of violent crimes from 2006 to 2015 in our data gives an observed value of 0.335, compared to a null mean of 0.0007 and standard error of 0.0127, which suggests a highly significant amount of spatial autocorrelation in violent crime totals.

In the next section, we develop several different Bayesian strategies for modeling violent crime over time and spatially between the areal neighborhoods of Philadelphia. We will fit our models on the violent crime data from 2006 to 2014, leaving data from 2015 for model comparison and evaluation.

3 Modeling Areal Crime Data over Space and Time

As described in Section 2, the areal units of our analysis are the 1336 US census block groups of Philadelphia (shown in Figure S1 in our supplementary materials (Balocchi and Jensen, 2019)).

For the remainder of this paper, we will use the terms “block group” and “neighborhood” interchangeably. The input data for our analysis is the number of violent crimes, $c_{it}$, reported in year $t$ within neighborhood $i$. Our temporal range is $t=1,\ldots,T$ with $T=10$, for the years 2006-2015 and our spatial range is $i=1,\ldots,n$ with $n=1336$, for all the block groups in Philadelphia.

As seen in the violent crime totals (averaged over time) in Figure S2 of our supplementary materials (Balocchi and Jensen, 2019), there are some substantial outlying neighborhoods with high violent crime totals relative to most of the city. These outliers (and general skewness in violent crime totals) motivates us to model violent crime totals on the logarithmic scale. This strategy has the additional benefit that linear changes over time in the logarithm of violent crime totals can be interpreted as percentage changes in raw violent crime totals.

However, because there are a small number of neighborhoods with zero crimes in some years, we need to consider a transformation that is defined at zero. Accordingly, we use the inverse hyperbolic sine transformation (Burbidge, Magee and Robb, 1988) that is centered to give values approximately equal to the logarithmic transformation. Specifically, we calculate our transformed violent crime totals as

[TABLE]

where $c_{it}$ is the total number of violent crimes reported in year $t$ within neighborhood $i$. A more common solution would be to add a small non-zero value to the counts, e.g. $\log(c_{it}+1)$. We prefer the inverse hyperbolic sine transformation as it is numerically equivalent to the $\log$ transformation for large counts but is a better approximation than the $\log(c_{it}+1)$ transformation for small counts.

An alternative modeling strategy for count data does not apply a transformation but assumes a Poisson distribution for the counts (Law, Quick and Chan, 2014; Li et al., 2014; Anderson and Ryan, 2017). The Poisson model would not work since our data is over-dispersed, and the more flexible negative-binomial distribution does not model mean and variance as intuitively as a normal model. In addition, the normal model is conjugate for the prior distributions we will be considering which eases posterior estimation.

3.1 Accounting for Neighborhood Level Covariates

We use a standard linear regression approach to account for the neighborhood-level economic, demographic and land use predictors of crime. Our transformed violent crime totals $y_{it}$ are modeled as,

[TABLE]

where $\mathbf{z}_{i}$ is the vector of predictor variables for neighborhood $i$ and $\boldsymbol{\gamma}$ is the vector of coefficients for those predictor variables, so $\mathbf{z}_{i}^{T}\boldsymbol{\gamma}=\sum_{d=1}^{6}\gamma_{d}z_{id}$.

As outlined in Section 2, we have $d=6$ predictor variables of crime for each neighborhood: population count, segregation, median household income, poverty, vacant proportion and commercial vs. residential proportion. We used square root transformations of vacant proportion, commercial vs. residential proportion and poverty and a logarithmic transformation of income to give a more linear relationship with the outcome variable.

Although yearly demographic and economic data is available after 2013, we avoid extrapolating values of the predictors to earlier years by modeling each predictor variable as static over the ten year period spanned by our crime data. We examine the estimated partial effects $\boldsymbol{\gamma}$ of these economic, demographic and land use predictors in Section 5.1.

Although there is interest in the partial effects of our crime predictors, our primary interest lies in the temporal trends captured by $e_{it}$ and the spatial correlation in these trends. With these time trends, we will be able to answer questions such as ‘what areas of the city are increasing or decreasing most quickly in terms of safety?’.

3.2 Time Trends with No Spatial Correlation

We can add a global linear trend over time into our model,

[TABLE]

where the scalar coefficient $\beta$ can be interpreted as the global percentage change in violent crime over time across the entire city of Philadelphia and $t$ takes on integer values from 1 to 10 to represent the years 2006-2015.

However, this model with only a global $\alpha$ and $\beta$ does not allow for heterogeneity between different neighborhoods in the overall level of violent crime or trend in violent crime over time. We can account for this heterogeneity through neighborhood-specific intercepts $\alpha_{i}$ and slopes $\beta_{i}$, which give us the model

[TABLE]

However, model (3.4) is over-parameterized: in fact, the effect of our static covariates is completely explained by the neighborhood-specific intercepts $\alpha_{i}$, so the same fit can be achieved by removing the covariates,

[TABLE]

Nonetheless, we can still estimate the partial effects of the covariates with an equivalent two-stage approach where we first fit $y_{it}=\alpha+\mathbf{z}_{i}^{T}\gamma+e_{it}$ and then fit the estimated residuals with the neighborhood-specific coefficient model, $\hat{e}_{it}=\alpha_{i}+\beta_{i}t+\epsilon_{it}$.

These neighborhood-specific model coefficients allow us to identify regions of Philadelphia with different levels of crime as well as different trends in crime over the past decade. This richer model is also motivated by fit to the data: a regression model with neighborhood-specific coefficients explains significantly more variation according to an F-test.

That said, we do not expect that every single neighborhood in Philadelphia would have unique coefficients, so we still risk over-parametrization with this model. We address this over-parameterization by imposing shared prior distributions for the neighborhood-specific coefficients from our time trend model (3.5),

[TABLE]

where we denote our collection of neighborhood specific coefficients with $\boldsymbol{\alpha}=(\alpha_{1},\ldots,\alpha_{n})$ and $\boldsymbol{\beta}=(\beta_{1},\ldots,\beta_{n})$. $\boldsymbol{\gamma}=(\gamma_{1},\ldots,\gamma_{d})$ collects the coefficients (partial effects) of the predictor variables which are shared by all neighborhoods.

We complete this model formulation by placing flat priors on the global means $\alpha_{0}$ and $\beta_{0}$, $p(\alpha_{0},\beta_{0})\propto 1$, and inverse gamma priors on the variance parameters

[TABLE]

The variance hyper-parameters are tuned in an empirical Bayes fashion so that the prior mean of the variance parameters is equal to the variance estimated from the model with no shrinkage, and the prior variance is small. Using non-informative priors for these variance parameters produced nearly identical results. See Section 4 of our Supplementary Materials (Balocchi and Jensen, 2019) for details.

This Bayesian hierarchical model shares information between neighborhoods by shrinking the neighborhood specific coefficients $\alpha_{i}$ and $\beta_{i}$ towards global parameters $(\alpha_{0},\beta_{0})$ for the entire city. For this reason, we refer to this approach as the global shrinkage model.

However, this global shrinkage model does not account for the spatial proximity between neighborhoods when sharing information. We expect close neighborhoods to behave similarly while we want distant neighborhoods to be informative but not as directly influential as adjacent ones. In other words, we may prefer a model that imposes local shrinkage rather than global shrinkage.

A model with local sharing of information would also be better able to address the substantial spatial correlation that we see in our application. Testing with Moran’s I shows that the residuals from the global shrinkage model are significantly spatially correlated. In the next subsection, we will explore conditional auto-regressive models for local sharing of information.

3.3 Time Trends with a Spatial Conditional Auto Regressive Model

A popular way of incorporating spatial information is through a prior distribution that is specified according to a Conditional Auto Regressive (CAR) model, which was introduced in its most general formulation by Besag (1974). The CAR model is a Gaussian Markov random field which induces spatial dependence through an adjacency matrix for the areal units, which in our case are neighborhoods in Philadelphia.

Several variations of this CAR framework are reviewed and compared in Lee (2011). In this paper, we will use the proper CAR formulation introduced by Leroux, Lei and Breslow (2000).

Let $\boldsymbol{\theta}$ denote a vector of elements that are potentially spatially correlated, such as our neighborhood-specific intercepts $\boldsymbol{\alpha}$ or slopes $\boldsymbol{\beta}$. Leroux, Lei and Breslow (2000) defines the distribution of each $\theta_{i}$ given the other $\theta_{-i}$ as a normal distribution centered at a weighted average of a global mean and the $\theta_{j}$’s from bordering neighborhoods,

[TABLE]

where $w_{ij}$ are adjacency weights that are equal to 1 if the neighborhoods $i$ and $j$ share a border and equal to 0 otherwise.

We collect these adjacency weights $w_{ij}$ into an adjacency matrix $\mathbf{W}$ that we assume (for now) to be known since we can easily use the shapefiles from the US Census Bureau to determine which of the 1336 neighborhoods (census block groups) share a border.

For now, we consider these adjacency weights $w_{ij}$ to be fixed. However, in Section 3.4 we will extend our model to allow those weights to vary since some borders may represent barriers between neighborhoods (e.g. highways or rivers), in which case we would not want to share information across that particular border.

The parameter $\rho\in[0,1]$ represents the strength of the spatial correlation between the components of $\boldsymbol{\theta}$, where larger values of $\rho$ correspond to a stronger influence of bordering neighborhoods. In the special case of $\rho=0$, the CAR prior (3.9) reduces to the global shrinkage prior (3.6)-(3.7).

It can be proved (Banerjee, Carlin and Gelfand, 2014, Ch.3) using Brook’s lemma (Brook, 1964), that the joint distribution of $\boldsymbol{\theta}$ is uniquely determined by the set of conditional distributions defined in 3.9:

[TABLE]

where $\mathbf{1}$ is a vector of 1’s and $\mathbf{D}_{W}-\mathbf{W}$ is the Laplacian matrix based on our neighborhood adjacency matrix $\mathbf{W}$. For values of $\rho$ in $[0,1)$ the joint distribution is proper, while for $\rho=1$ the distribution is degenerate (Lee, 2011). By adding the constraint $\sum_{i}(\theta_{i}-\theta_{0})=0$ we can get a distribution for a $n$-dimensional vector, concentrated in a $(n-1)$-dimensional subspace; this is known as the intrinsic CAR by Besag, York and Mollié (1991).

We will employ this CAR model as prior distributions for the vectors of time trend coefficients $\boldsymbol{\alpha}$ and $\boldsymbol{\beta}$. We assume $\boldsymbol{\alpha}$ and $\boldsymbol{\beta}$ are a priori independent. In vector form, the CAR model (3.9) corresponds to the following prior distributions for $\boldsymbol{\alpha}$ and $\boldsymbol{\beta}$,

[TABLE]

where $\boldsymbol{\Sigma}^{-1}=\rho(\mathbf{D}_{W}-\mathbf{W})+(1-\rho)\mathbf{I}$.

We use the same prior distributions for $\alpha_{0}$ and $\beta_{0}$ and our variance parameters as in the global shrinkage model in the previous subsection. For the additional spatial parameter $\rho$, we choose a Beta$(10,10)$ prior distribution which has mean equal to $0.5$ and a small variance in order to avoid the endpoints of the interval $[0,1]$.

The posterior distributions for the spatial CAR model and the global shrinkage model (Section 3.2) can be implemented via a Gibbs sampler (Geman and Geman, 1984). Implementation details are given in our supplementary materials (Balocchi and Jensen, 2019).

3.4 Allowing Neighborhood Border Weights to Vary

For most types of areal data, the weights $\mathbf{W}$ that encode the spatial connection between the areal units are considered to be fixed and known. In our data context, the areal units are neighborhoods and the weights $\mathbf{W}$ encode which neighborhoods share a border and hence induce shrinkage on each other in our spatial CAR models outlined in Section 3.3.

However, within any large city, some borders between neighborhoods consist of natural or artificial barriers such as rivers, highways or train tracks. These barriers could reduce the similarity in crime trends between neighborhoods, and so we would not want to shrink estimates across those barriers. The implication of these barriers for the spatial CAR models in Section 3.3 are that some weights $w_{ij}=1$ should really be $w_{ij}=0$ since those neighborhoods share a border that is actually a barrier.

Attempting to set which borders should actually be barriers manually would be tedious for a large city and also require extensive domain knowledge and subjective decision making. We instead prefer to infer these barriers from the data by allowing a subset of weights $w_{ij}$ to be random variables in our model.

Specifically, we consider the set of indices of pairs of neighborhoods which share a border according to the geography of Philadelphia. The matrix $\mathbf{W}$ is symmetric so the random variables $w_{ij}$ and $w_{ji}$ are considered to be the same object. We model the $w_{ij}$ for neighborhood pairs that share a border as Bernoulli random variables with an prior probability $\phi$ of $w_{ij}=1$. Any weights $w_{ij}=0$ according to the geography of Philadelphia will remain fixed at $w_{ij}=0$ since we do not want to form connections between non-proximal neighborhoods.

We expect a priori that the probability $\phi$ will be close to 1, since relatively few borders between neighborhoods actually should be barriers. For this reason we choose the prior for $\phi$ to be a ${\rm Beta}(9,1)$ distribution which has mean close to one and small variance.

Moreover, we expect that the spatial distribution of the neighborhood-specific crime levels $(\alpha_{i})$ may be different from the neighborhood-specific trends in crime over time $(\beta_{i})$, so we allow for different barriers when we model the distribution of $\boldsymbol{\alpha}$ and $\boldsymbol{\beta}$. In particular, we consider two random matrices $\mathbf{W}^{\alpha}$ and $\mathbf{W}^{\beta}$ where a subset of the elements of these matrices are random as described above: $w_{ij}^{\alpha}\,|\,\phi^{\alpha}\sim{\rm Bernoulli}(\phi^{\alpha})$ and $w_{ij}^{\beta}\,|\,\phi^{\beta}\sim{\rm Bernoulli}(\phi^{\beta})$ for neighborhood pairs $(i,j)$ that share a border.

These two weight matrices then determine the local shrinkage of our spatial CAR model from the previous subsection:

[TABLE]

where $\boldsymbol{\Sigma}_{\alpha}^{-1}=\rho\cdot(\mathbf{D}_{W^{\alpha}}-\mathbf{W}^{\alpha})+(1-\rho)\mathbf{I}$ and $\boldsymbol{\Sigma}_{\beta}^{-1}=\rho\cdot(\mathbf{D}_{W^{\beta}}-\mathbf{W}^{\beta})+(1-\rho)\mathbf{I}$.

Allowing variable border weights can lead to over-parametrization since we are adding as many parameters as the number of borders, which makes the shrinkage imposed by prior parameters $\phi^{\alpha}$ and $\phi^{\beta}$ important. A more sophisticated approach, which is the focus of ongoing work, would be to partition our areal units into clusters with barriers represented as cluster boundaries.

To implement this extended model with some variable border weights, a step is added to our Gibbs sampler that samples each border weight conditional on the current values of the other model parameters. Details are given in our supplementary materials (Balocchi and Jensen, 2019).

The idea of detecting discontinuities at boundaries is often referred to as wombling after the seminal work of Womble (1951) and has been very popular in the disease mapping literature. However, most papers have approached detection of boundaries as a selection problem that is performed after inference (see, e.g. Boots, 2001; Li, Banerjee and McBean, 2011; Banerjee et al., 2012; Lu and Carlin, 2005; Lee and Mitchell, 2013).

In contrast, we incorporate the possibility of discontinuities at boundaries directly into our model through variable $\mathbf{W}^{\alpha}$ and $\mathbf{W}^{\beta}$, which allows us to incorporate potential barriers into our estimation of neighborhood-specific parameters $\boldsymbol{\alpha}$ and $\boldsymbol{\beta}$. Lee and Mitchell (2012) and Lu et al. (2007) take a similar approach in the context of disease mapping, but with a more elaborate model for ${\rm P}(w_{ij}=1)$ that is a function of dissimilarity between covariate values in units $i$ and $j$.

In Section 4, the different models presented in this section are compared in terms of their accuracy of their in-sample and out-of-sample predictive accuracy. We then visualize the estimated trends in crime over time in Philadelphia and discuss several insights from our results in Section 5.

4 Comparison of Predictive Accuracy

In the previous section, we outlined a no shrinkage model (Section 3.2) and several hierarchical Bayesian models for estimating the neighborhood-level trend in crime over time, including a global shrinkage model (Section 3.2), a spatial CAR models for local shrinkage (Section 3.3), and finally an extension of the spatial CAR model to allow a subset of border weights to vary (Section 3.4).

We now compare each of these model alternatives based on several measures of the accuracy of their predictions on both in-sample and out-of-sample hold-out data. Recall that we have 10 years of crime data for the city of Philadelphia, from the beginning of 2006 to the end of 2015. We estimate each model using the crime data for the first nine years (2006-2014).

We assess the in-sample accuracy of each model by computing the mean squared error of the predictions of violent crime totals for 2014, which is a year that was included in model estimation,

[TABLE]

We assess the out-of-sample accuracy of each model by computing the mean squared error of the predictions of violent crime totals for 2015, which is a year that was not included in model estimation,

[TABLE]

To ensure our evaluation is not overly dependent on any idiosyncratic aspects of the 2015 data, we also calculate the cross-validated out-of-sample accuracy of each model by calculating the mean square error ${\rm MSE}_{\rm out}^{t}$ when using year $t$ as the hold out data in the same way that 2015 is used as the hold out data in (4.2), i.e.

[TABLE]

In Table 1, we compare the predictive accuracy of four different models with neighborhood-specific coefficients outlined in Section 3: 1. the time trend model (3.5) without shrinkage between neighborhoods, 2. the global shrinkage model with priors (3.6) and (3.7), 3. the local shrinkage model with spatial CAR priors (3.11) and (3.12) and 4. the local shrinkage spatial CAR model with variable borders (3.13) and (3.14). For additional reference, we also provide the mean square error for fitting a single trend (“Global $\alpha\,,\,\beta$”) across the entire city.

We see in Table 1 that the model with a global trend over time (“Global $\alpha\,,\,\beta$”) for the entire city has very poor predictive accuracy compared to the models that allow neighborhood-specific time trends (“Separate $\alpha_{i}\,,\,\beta_{i}$”).

Among the neighborhood-specific time trend models, the global shrinkage model has substantially lower out-of-sample mean square errors than the baseline time trend model without any shrinkage between neighborhoods. The best in-sample mean squared error was achieved by the model without shrinkage, as we expect from the least square method, though at a cost of having the worst out-of-sample accuracy.

The model with local shrinkage via the spatial CAR prior further reduces the out-of-sample mean square errors compared to the global shrinkage model. The model that allows variable borders does not further improve the out-of-sample mean squared errors, though we explore in Section 5.3 that it helps with the interpretation.

Table 1 also provides Moran’s $I$ measure of spatial autocorrelation, calculated on the posterior mean of the neighborhood-specific time trends ($\beta_{i}$’s). We see that the spatial CAR model induces a larger spatial correlation in the $\beta_{i}$’s than the models with global shrinkage or without shrinkage. The local shrinkage model has a Moran’s I value of 0.61 (s.e. = 0.016), which suggests there is substantial spatial autocorrelation in the change in crime within Philadelphia.

In summary, allowing for local shrinkage of the neighborhood-specific crime trend coefficients via the spatial CAR priors (3.11) and (3.12) leads to the best out-of-sample predictive accuracy. In Section 5, we visualize the parameters of this model and discuss the implications of these results for crime in Philadelphia.

Although the variable border model extension does not improve out-of-sample predictive accuracy, we will also see in Section 5 that visualizing the borders that have been turned into barriers by this model provide insight into discontinuities in crime trends in the city of Philadelphia.

5 Interpretation of Model Parameters

In Section 5.1, we examine the estimated partial effects for the static predictor variables created from the data outlined in Section 2. We then visualize and compare the estimated neighborhood-specific levels ($\alpha_{i}$’s) and time trends ($\beta_{i}$’s) on crime from our different models in Section 5.2. In Section 5.3, we examine the results from our model extension outlined in Section 3.4 that allows a subset of neighborhood borders in Philadelphia to be estimated as barriers. Finally, in Section 5.4 we discuss the neighborhoods with the most extreme levels and changes in crime over time over the past ten years in Philadelphia.

5.1 Partial Effects of Static Predictors

Figure 2 gives the estimated partial effects $\gamma_{d}$ for each static predictor variable $d$ from the four models outlined in Section 3. We provide additional numerical details in Table S1 of our supplementary materials (Balocchi and Jensen, 2019).

We see that among the six predictor variables created in Section 2, only the segregation measure is not a significant predictor of crime. All predictor variables are on the same scale and so we can directly compare the values of their partial effects.

We see that the strongest predictors of crime are total population and the commercial versus residential proportion, with more populated and more commercial neighborhoods being associated with higher crime. Income and poverty are also significantly predictive of violent crimes but we must be more cautious about interpreting these partial effects given the high collinearity between income and poverty. Each of these observations on the partial effects $\boldsymbol{\gamma}$ is relatively consistent across the four models outlined in Section 3.

5.2 Visualizing Neighborhood-Specific Coefficients

Our primary interest in terms of interpretation are the estimated neighborhood-specific coefficients, $\alpha_{i}$’s and $\beta_{i}$’s, that represent the level of violent crimes and change in violent crimes over time in Philadelphia, respectively.

In Figures 3 and 4, we give maps where each block group in Philadelphia is colored by the estimated neighborhood-specific levels of crime $\hat{\alpha}_{i}$ and changes in crime over time $\hat{\beta}_{i}$ respectively, from the four models outlined in Section 3. We see substantial heterogeneity between neighborhoods in Philadelphia, both in terms of the their estimated crime levels ($\hat{\alpha}_{i}$’s) and changes in crime over time ($\hat{\beta}_{i}$’s). Regardless of the model, most neighborhoods in the city show decreasing trends in crime over time (negative $\beta$’s) with a small subset of neighborhoods showing an increasing trend.

The shrinkage imposed by the global shrinkage model is more visually striking for the change in violent crime over time than the overall level of crime. The maps of the $\hat{\alpha}_{i}$’s from the no shrinkage and global shrinkage models are almost indistinguishable in Figure 3 whereas the map of the $\hat{\beta}_{i}$’s from the global shrinkage model has been shifted substantially compared to the no shrinkage map in Figure 4. This observation suggests that there is more substantial heterogeneity between neighborhoods in terms of their overall level of crime compared to their change in crime over time.

This heterogeneity in the mean level of crime is expected as it is influenced by many years of transformation in the city of Philadelphia that led to its current built and social environment. Differences in these overall spatial crime patterns can be addressed by urban planners, whose effects are long-lasting (Johnson et al., 2008). In contrast, differences in the trend over time identify shorter-term patterns, which can be addressed with interventions by local police departments.

The overall level of crime also seems to have a greater inherent spatial correlation between proximal neighborhoods than the change in crime over time. The Moran’s I values calculated from the estimated $\hat{\alpha}_{i}$’s are I = 0.33 for both the no shrinkage and global shrinkage models, compared to the value of I = 0.17 from the estimated $\hat{\beta}_{i}$’s for those same models in Table 1. This is clear also from the maps from the no shrinkage model (top left) in Figures 3 and 4: the estimated $\hat{\beta}_{i}$’s are more “spotty” and less smooth than the corresponding map of the $\hat{\alpha}_{i}$’s.

However, once we build spatial correlation into our model via the spatial CAR prior (3.11)-(3.12), the resulting $\hat{\beta}_{i}$’s are more spatially correlated than the resulting $\hat{\alpha}_{i}$’s, as can be seen in the lower left of Figures 3 and Figures 4 as well as the corresponding Moran’s I = 0.53 for the $\hat{\alpha}_{i}$’s versus I = 0.61 for the $\hat{\beta}_{i}$’s. Note that all these reported Moran’s I values have a standard error approximately equal to 0.016, and so they are all significantly different from the null hypothesis of no spatial autocorrelation.

Although the smoother maps from the spatial CAR model (lower left of Figures 3 and  4) ease interpretation by identifying larger regions of the city with similar crime dynamics, there is the potential to over-shrink certain neighborhoods that should actually stand out from their neighbors. In any large city, natural or artificial barriers such as rivers, highways or rail lines create discontinuities between neighborhoods which should not be smoothed over. In Section 5.3, we examine the results from our model extension that allows a subset of borders between neighborhoods to be turned into barriers.

5.3 Borders turned into Barriers

In Section 3.4, we extended the spatial CAR model to allow a subset of the weights $w_{ij}$ to vary, which allows the borders ($w_{ij}=1$) between some neighborhoods to be changed into barriers ($w_{ij}=0$); the latter prevent shrinkage between two bordering neighborhoods. Our model has separate weight matrices $\mathbf{W}^{\alpha}$ and $\mathbf{W}^{\beta}$, so a particular border can be turned into a barrier either for the level of crime ($\alpha_{i}$’s) or the change in crime over time ($\beta_{i}$’s) or both. Using this model, we estimate the posterior probability that we change a border into a barrier for each border between proximal neighborhoods in Philadelphia.

Figure 5 gives the distribution of the estimated posterior probability of a border being turned into barrier for each border encoded in the weight matrices $\mathbf{W}^{\alpha}$ and $\mathbf{W}^{\beta}$. These distributions seem to have two components: a main mode representing the behavior of the majority of the borders, which has a low probability of being turned into a barrier, and a “tail” component which has a higher probability of being turned into a border.

It is clear that many more borders have a high probability of being a barrier for the level of crime ($\alpha_{i}$’s) compared to the change in crime over time ($\beta_{i}$’s). In other words, our variable border model is detecting more discontinuities between bordering neighborhoods in the level of crime compared to the change in crime over time. In Section 5 of our Supplementary Materials (Balocchi and Jensen, 2019), we explore an alternative model that only allows variable borders for the mean level of crime.

In the lower right panels of Figures 3 and  4, we provide maps of Philadelphia where we have highlighted any borders between neighborhoods that have been inferred by our model to have a high probability of being barriers. These particular highlights are based on posterior probabilities larger than $60\%$ for $\mathbf{W}^{\alpha}$ and larger than 50% for $\mathbf{W}^{\beta}$.

We see in the lower right panel of Figure 3 that barriers have been detected around several parks including Fairmount Park, Wissahickon Valley Park, and Pennypack Creek Park (indicated by the black numbers 1, 2 and 3 respectively in the lower right panel of Figure 3). In these cases, our model has automatically detected several natural geographic structures within Philadelphia as locations which have discontinuities in the level of crime.

We also see that some estimated barriers have isolated particular neighborhoods from their proximal neighbors. For example, the neighborhood of Bridesburg (indicated by the black number 4 in the lower right panel of Figure 4) seems to have a much more positive trend on crime over time than its surrounding neighborhoods.

As barriers highlight the boundaries of regions that display differences in either in the level of crime or the trend in crime over time, these barriers can be used by police departments and city planners for delineating the possible limits of effectiveness for interventions or as potential targets for interventions themselves.

5.4 Neighborhoods with Most Extreme Crime Trends

To further understand which regions of Philadelphia have the most extreme levels of crime and trends in crime over time, we can examine the most extreme intercepts ($\alpha_{i}$’s) and slopes ($\beta_{i}$’s) found by our fitted models. Specifically, we focus on the estimated $\alpha_{i}$’s and $\beta_{i}$’s from the local shrinkage spatial CAR model (3.11)-(3.12) that had the best out-of-sample predictive performance in Table 1.

Figure 6 provides maps that highlight the most extreme (largest 50 and smallest 50) neighborhoods in terms of the estimated level of crime ($\hat{\alpha}_{i}$’s) and in terms of the estimated change in crime over time ($\hat{\beta}_{i}$’s).

We see that the region of University City in West Philadelphia (black number 1 in the top panel of Figure 6) is an interesting transitional area that contains both neighborhoods with the highest and lowest levels of crime in the city. We also see that the area of Frankford (black number 2 in the top panel of Figure 6) has neighborhoods with high levels of crime. This area is a major transportation hub for the Northeast region of Philadelphia.

The SW region of Philadelphia, specifically the Elmwood and Eastwick neighborhoods (black number 3 in the bottom panel of Figure 6) have seen some of the largest reductions in crime over the past decade in Philadelphia. We also see some regions of the city that are showing increases in crime over that same time period, such as the Wissinoming and Tacony neighborhoods (black number 4 in the bottom panel of Figure 6) that are just to the northeast of the high crime neighborhoods of Frankford (black number 2 in the top panel of Figure 6).

In Section 5 of our Supplementary materials (Balocchi and Jensen, 2019), we provide additional visualizations of the neighborhood-specific parameters that are significantly different from the overall mean across the city as well as the widths of the credible intervals for these parameters.

6 Discussion

Reliable estimation of the change in crime over time at the local neighborhood level is a crucial step towards a better understanding of the determinants of public safety in large urban areas. With a focus on the city of Philadelphia, we have explored several Bayesian approaches to modeling crime trends within the areal units of neighborhoods while sharing information either globally or locally across the city.

Imposing local shrinkage between proximal neighborhoods via a spatial conditional autoregressive (CAR) prior gives the best out-of-sample predictions of violent crime compared to models that impose global shrinkage or no shrinkage at all between neighborhoods. We also explore allowing the weight matrix of our spatial CAR model to vary in order to detect neighborhood borders that represent spatial discontinuities in the level of crime or change in crime over time. In this way, we automatically detect several natural barriers in the geography of Philadelphia. Our model estimates also identify the regions of Philadelphia with the most extreme levels of violent crime as well as the largest increases and reductions in crime over the period of 2006-2015.

7 Acknowledgements

We thank Rachel Thurston and Theresa Smith for helpful contributions to our modeling efforts and interpretation of our results.

{supplement}

[id=suppA] \stitleSupplementary materials for “Spatial modeling of trends in crime over in Philadelphia” \slink[doi]COMPLETED BY THE TYPESETTER \sdatatype.pdf \sdescriptionWe provide maps outlining the block group structure, population count, and distribution of violent crimes in Philadelphia. We give implementation details for our fitted models. We compare results under an alternative choice of prior distributions and under an alternative model that only allows variable borders for the mean level of crime. We provide additional numerical details about our estimated partial effects. We also provide visualizations of the significance and widths of credible intervals for the neighborhood-specific parameters.

8 Maps of Data in Philadelphia

Figure S1 (left) gives a map outlining the 1336 block groups in Philadelphia. Figure S1 (right) shows population count for each block group in Philadelphia.

In Figure S2 (left), we give the count of violent crimes per year in each block group, averaged over the years 2006-2015. We see substantial heterogeneity across block groups in the average counts of violent crimes per year. There are several outlying values: particular block groups that have much higher average violent crime counts. The largest among these is the Market East neighborhood in central Philadelphia.

These outlying neighborhoods motivate us to examine violent crime totals on the log scale. In Figure S2 (right), we give the average of the logarithm of the count of violent crimes per year in each block group, averaged over the years 2006-2015. We can see more details of the spatial distribution of violent crime on the log scale. Modeling crime on the log scale has the additional benefit that changes in log crime can be interpreted as percentage changes in crime.

In both Figure S2 (left) and (right), we see evidence of spatial correlation in violent crime totals between proximal block groups throughout the city.

9 Gibbs sampling

In section 3 of our paper, we described the different models considered in this work; here we illustrate the Gibbs sampling strategy to sample from their posterior distributions. We are mainly interested in the coefficients $\boldsymbol{\gamma}$ for the predictors and in the collection of neighborhood-specific coefficients $(\boldsymbol{\alpha},\boldsymbol{\beta})$. $\boldsymbol{\gamma}$ have a Normal prior distribution with mean zero and covariance matrix proportional to the identity matrix. $\boldsymbol{\alpha}$ and $\boldsymbol{\beta}$ also have Normal prior distributions with mean zero, but have different covariance matrices depending on the model. In the global shrinkage model, the covariance matrices are proportional to the identity matrix, while in the local shrinkage CAR model the covariance matrices depend on the Laplacian matrix of the geography.

We denote with $\mathbf{Y}$ the $N$-dimensional vector (where $N=nT$) found by concatenating the $y_{it}$’s, ordered by block group: $\mathbf{Y}=(y_{11},\ldots,y_{1T},\ldots,y_{n1},$ $\ldots,y_{nT})$; moreover let $\boldsymbol{\theta}=(\boldsymbol{\gamma},\boldsymbol{\alpha},\boldsymbol{\beta})$ represent the collection of the coefficients. Let the matrix $\mathbf{X}$ be designed so that the covariates corresponding to block group $i$ at time $t$ are contained in row $(i-1)T+t$ and by multiplying this row with the vector of parameters we get $X_{(i-1)T+t}^{\gamma}\boldsymbol{\theta}=\sum_{j}\gamma_{j}z_{ij}+\alpha_{i}+\beta_{i}t$. With this notation, the conditional distribution of the data is given by $\mathbf{Y}|\boldsymbol{\theta},\sigma^{2}\sim N(\mathbf{X}\boldsymbol{\theta},\sigma^{2}\mathbf{I})$.

We can also combine the prior distributions of $\boldsymbol{\gamma},\boldsymbol{\alpha},\boldsymbol{\beta}$ to get the distribution of $\boldsymbol{\theta}$: let $\boldsymbol{\theta}_{0}=(\mathbf{0},\alpha_{0}\mathbf{1},\beta_{0}\mathbf{1})$ a ($d+2n$)-dimensional vector representing the conditional mean of $\boldsymbol{\theta}$ and let $\boldsymbol{\Omega}_{0}=\boldsymbol{\Sigma}_{0}^{-1}$ be the block matrix representing its precision matrix. Since $p(\boldsymbol{\gamma})={\rm N}(\mathbf{0},\tau_{\gamma}^{2}\cdot\mathbf{I})$, the first $d\times d$ diagonal block of $\boldsymbol{\Omega}_{0}$ is equal to $\tau_{\gamma}^{-2}\mathbf{I}$; the next two $n\times n$ diagonal blocks instead are the precision matrices of $\boldsymbol{\alpha}$ and $\boldsymbol{\beta}$: $\tau_{\alpha}^{-2}\boldsymbol{\Sigma}^{-1}$ and $\tau_{\beta}^{-2}\boldsymbol{\Sigma}^{-1}$, where $\boldsymbol{\Sigma}=\mathbf{I}$ in the global shrinkage model and $\boldsymbol{\Sigma}^{-1}=[\rho(\mathbf{D}_{W}-\mathbf{W})+(1-\rho)\mathbf{I}]$ in the spatial CAR model; the remaing blocks are zero matrices. Then $p(\boldsymbol{\theta}|\boldsymbol{\theta}_{0},\tau_{\alpha}^{2},\tau_{\beta}^{2})\propto\exp\left(-\frac{1}{2}(\boldsymbol{\theta}-\boldsymbol{\theta}_{0})^{T}\boldsymbol{\Omega}_{0}(\boldsymbol{\theta}-\boldsymbol{\theta}_{0})\right)$. Moreover, we set non-informative flat priors on $\alpha_{0},\beta_{0}$, and the variance hyper-parameters $\sigma^{2}$, $\tau_{\gamma}^{2}$, $\tau_{\alpha}^{2}$ and $\tau_{\beta}^{2}$ have Inverse-Gamma priors, tuned in an Empirical Bayes fashion.

Posterior conditional distribution of $\boldsymbol{\theta}$

With this notation, we can find the conditional posterior distribution of $\boldsymbol{\theta}$ as we would do in a usual linear regression:

[TABLE]

where

[TABLE]

Posterior conditional distribution of $\boldsymbol{\theta}_{0}$

Similarly, the posterior distribution of the mean hyper-parameters $\alpha_{0}$ and $\beta_{0}$ can be found as

[TABLE]

Posterior conditional distribution of $\sigma^{2}$, $\tau^{2}_{\gamma}$, $\tau^{2}_{\alpha}$ and $\tau^{2}_{\beta}$

For the variance hyper-parameters $\sigma^{2},\tau^{2}_{\gamma},\tau_{\alpha}^{2}$ and $\tau_{\beta}^{2}$, the prior distributions are

[TABLE]

where the hyper-parameters are tuned in an empirical Bayes fashion so that the prior mean of the variance parameters is equal to the variance estimated from the model with no shrinkage, and the prior variance is small. The conditional posterior distributions are also Inverse-Gamma:

[TABLE]

Posterior conditional distribution of $\rho$

The prior distribution on $\rho$ is Beta$(10,$ $10)$, and since its conditional posterior distribution does not have a closed form, we sample this parameter with a Metropolis Hasting procedure. Given its past value $\rho^{t-1}$ we propose a new candidate $\rho^{*}$ with density $g(\rho^{*}|\rho^{t-1})={\rm Beta}(b\rho^{t-1}/(1-\rho^{t-1}),b)$; this parametrization allows the mean to be $\rho^{t-1}$ and the variance to be small when we choose $b=10$. The acceptance probability is then

[TABLE]

where the posterior conditional distribution $p(\rho|e.e.)$ is proportional to the product of the prior distribution of $(\boldsymbol{\alpha},\boldsymbol{\beta})$ given $\rho$ and the prior of $\rho$. We use the notation $e.e.$ to denote “everything else”, i.e. the current values of all other parameters in the model.

Posterior conditional distribution of $\mathbf{W}$

Finally, in model (3.13)-(3.14) we allow the adjacency matrix itself to be random. We consider all the pairs of regions that share a border ($w_{ij}=1$) and we allow those borders to potentially become barriers ($w_{ij}=0$). We model these variable weights as $w_{ij}^{\alpha}|\phi^{\alpha}\sim{\rm Bern}(\phi^{\alpha})$ and independently, $w_{ij}^{\beta}|\phi^{\beta}\sim{\rm Bern}(\phi^{\beta})$, with $\phi^{\alpha},\phi^{\beta}\overset{iid}{\sim}{\rm Beta}(1,9)$.

With these prior distributions, the conditional posterior distribution for $\mathbf{W}^{\alpha}$ is

[TABLE]

Note that, because of the determinant term, the entries of $\mathbf{W}^{\alpha}$ are not independent a posteriori. Thus we sample each entry $w_{ij}^{\alpha}=w_{ji}^{\alpha}$ conditional on the rest of the matrix $\mathbf{W}^{\alpha}_{-ij}$ as $p(w_{ij}^{\alpha}=1|e.e.)=q$, where

[TABLE]

A highly similar procedure (with the obvious substitutions) is used to sample the entries of $\mathbf{W}^{\beta}$.

Posterior conditional distribution of $\phi$

To express the prior information that only a small percentage of the borders should be turned into barrier, the prior distribution of $\phi^{\alpha}$ and $\phi^{\beta}$ is Beta$(1,9)$. Since the $w_{ij}^{\alpha}$ and $w_{ij}^{\beta}$ are Bernoulli distributed, the posterior distributions for $\phi^{\alpha}$ and $\phi^{\beta}$ are

[TABLE]

10 MCMC Implementation Details

The results reported in Table 1 of our paper for the hierarchical models have been implemented using Gibbs sampling. In particular, for each model 1000 samples where used, after discarding a burn-in period of 50 iterations and thinning every 2 samples. By running multiple chains and superimposing their trace plots, we noted that the convergence happened after a relatively short time and that samples were not highly correlated.

11 Prior Robustness for Variance Hyperparameters

In Section 3 of our main paper, we used priors for the variance parameters with hyper-parameters that were tuned in an Empirical Bayes fashion. In this section, we show that highly similar results are obtained when using more non-informative prior distributions on these variance parameters.

In particular, we consider a uniform prior on $\boldsymbol{\gamma},\log\sigma,\tau_{\alpha}$ and $\tau_{\beta}$, which is equivalent to:

[TABLE]

Table S1 is equivalent to Table 1 in our main paper but with results from the estimated models that use the non-informative priors given above. Only the global, spatial CAR and variable border model results are reported since priors are not involved in the no-shrinkage model (3.5).

Examining Table S1, we see almost the exact same predictive results as the predictive results given in Table 1 of our main paper. There are very slight numerical differences in the Spatial CAR in-sample error, the Global Shrinkage out-of-sample error and the Variable Borders Moran’s I, but these differences could easily be attributed to MCMC sampling variability.

12 Additional Model Results

In Section 12.1 we report the numerical estimates of the partial effects, which are shown in Figure 2 of the main paper. In Section 12.2 we describe the results from a model with variable borders for $\alpha_{i}$ but fixed borders for $\beta_{i}$. In section 12.3 we provide different visualizations of the variability in the estimated neighborhood-specific coefficients $\alpha_{i}$ and $\beta_{i}$.

12.1 Numerical results for partial effects

In table S2 we report the numerical values corresponding to the partial effects shown in Figure 2 of the main paper. For each model with neighborhood-specific coefficients outlined in Section 3 of the paper, we report maximum likelihood estimates, standard errors, posterior means and posterior standard deviations.

12.2 Random borders for only $\boldsymbol{\alpha}$

Given the smaller number of barriers detected for the $\beta_{i}$’s compared to the $\alpha_{i}$’s in Figure 5 of our main paper, we also implemented an alternative model where the adjacency structure $\mathbf{W}^{\beta}$ for the $\beta_{i}$’s is considered fixed, and only the adjacency matrix $\mathbf{W}^{\alpha}$ for the $\alpha_{i}$’s is allowed to vary.

The in-sample MSE of 0.0711 for the model with variable $\mathbf{W}^{\alpha}$ and fixed $\mathbf{W}^{\beta}$ is worse then the in-sample MSE of 0.0706 for the model with variable $\mathbf{W}^{\beta}$ and $\mathbf{W}^{\alpha}$. However, the variable $\mathbf{W}^{\alpha}$ and fixed $\mathbf{W}^{\beta}$ model does have a slightly better out-of-sample MSE of 0.1050 compared to the out-of-sample MSE of 0.1069 for the model with variable $\mathbf{W}^{\beta}$ and $\mathbf{W}^{\alpha}$.

These results provide a further indication that there is stronger signal in the data for detecting discontinuities for the mean level of crime between neighborhoods compared to discontinuities for the trend in crime over time between neighborhoods. However, we still report the results for the model with variable $\mathbf{W}^{\alpha}$ and $\mathbf{W}^{\beta}$ in our main paper, as it provides additional insight and interpretation in Figure 4 of the main paper.

12.3 Variability of neighborhood-specific coefficients

In Section 5.4 of our paper, we report the neighborhoods with the highest and lowest estimates of the mean level of crime $\alpha_{i}$ and time trend $\beta_{i}$. As a supplement to these results, in Figure S3 we map the Philadelphia neighborhoods which are “significant”, in the sense that their $95\%$ credible intervals do not contain the global mean $\bar{\alpha}$ or $\bar{\beta}$ across all neighborhoods in the city.

In Figure S3, we see many more neighborhoods with significant differences in terms of their $\alpha_{i}$’s, which is another indication that the variation in the mean level of crime is larger than the variation in the time trend in crime ($\beta_{i}$’s). In the plot for $\beta_{i}$’s, we find a smaller number of neighborhoods with values that are significantly different than the overall mean, but the existence of these neighborhoods confirms the presence of the space-time interaction found in previous studies (Law, Quick and Chan, 2014; Li et al., 2014).

In Figure S4, we visualize the width of the $95\%$ credible intervals for each neighborhood-specific $\alpha_{i}$ and $\beta_{i}$. It is interesting to observe that the interval widths are smallest for areal units that border many other units and largest for areal units with very few neighbors. The neighborhoods with smallest widths are the parks (Fairmount, Wissahickon and Pennipack) which are surrounded by many block groups due to their large surface area. The neighborhoods with largest widths are at the border of the city with only one or two neighboring units. This phenomenon is more striking for the trends in crime over time ($\beta_{i}$’s) for which there is a less strong signal in the data compared to the mean level of crime ($\alpha_{i}$’s) for each neighborhood.