Joint Modeling of Multiple Crimes : A Bayesian Spatial Approach

A multivariate Bayesian spatial modeling approach was used to jointly model the counts of two types of crime, i.e., burglary and non-motor vehicle theft, and explore the geographic pattern of crime risks and relevant risk factors. In contrast to the univariate model, which assumes independence across outcomes, the multivariate approach takes into account potential correlations between crimes. Six independent variables are included in the model as potential risk factors. In order to fully present this method, both the multivariate model and its univariate counterpart are examined. We fitted the two models to the data and assessed them using the deviance information criterion. A comparison of the results from the two models indicates that the multivariate model was superior to the univariate model. Our results show that population density and bar density are clearly associated with both burglary and non-motor vehicle theft risks and indicate a close relationship between these two types of crime. The posterior means and 2.5% percentile of type-specific crime risks estimated by the multivariate model were mapped to uncover the geographic patterns. The implications, limitations and future work of the study are discussed in the concluding section.


Introduction
Ecological studies of crime are of great interest to geographers and criminologists as they can reveal the geographic pattern of crime risks as well as the relevant risk factors explaining that pattern.Researchers have used non-spatial regression models in crime analysis.For example, Craglia et al. [1] adopted a standard logistic regression model to explore the relationship between socioeconomic characteristics and high-intensity crime areas.Such non-spatial models assume that the observations are independent and identically distributed.This assumption, however, is usually problematic in empirical crime studies.In practice, adjacent areas often have similar crime data, i.e., crime data are very likely to be spatially auto-correlated, especially at small-area scales as often indicated by clusters appearing on crime maps.Ignoring this spatial dependence in ecological studies can result in estimates with underestimated standard errors [2].Thus, methodologically, the traditional non-spatial models are inappropriate to analyze crime at the small-area level.
Spatial autocorrelation is typically modeled using the spatial lag model [3], the spatial error model [4] or other frequentist statistical models [5] in crime analysis.However, there are limitations with frequentist models.The spatial lag and error models [6] can only be used to model continuous variables (i.e., outcomes of rates, not counts) following a normal distribution.In addition, frequentist spatial models cannot process available information systematically [7].These approaches are also ISPRS Int.J. Geo-Inf.2017, 6, 16 2 of 16 incapable of dealing with the small number problem, which may produce unreliable risk estimations if not accounted for in the model [8].
The limitations of the frequentist models in small-area crime analysis can be overcome by Bayesian spatial models.Bayesian statistics, unlike frequentist statistics, treats unknown parameters as random variables expressed in terms of probabilities.It combines prior knowledge and observation data to obtain posterior distributions for parameters of interest.The implementation of Markov chain Monte Carlo (MCMC) methods and software such as WinBUGS to implement MCMC [7,9] make it possible to fit complex spatial models.Although Bayesian spatial models can tackle problems encountered in frequentist approaches, and thus are advantageous in small-area analysis, they are still seldom applied in crime studies.This is changing, however, as more researchers have become aware of the advantages of Bayesian spatial models.A small but growing body of research has adopted Bayesian spatial modeling approaches in crime analysis [10][11][12][13][14][15][16][17].For example, researchers using Bayesian spatial models have investigated contextual influences on domestic violence [12,15,18], juvenile offenders [14] and property crime [13].
The existing criminological research conducted by using Bayesian spatial models typically focuses on one type of crime.In these applications, the modeling methods adopted are univariate in nature.Counts of different types of crime can be interrelated and are multivariate to some extent as the crimes may share some missing or unmeasured variables.Modeling interrelated crime counts by aggregating them all together or separately using a univariate approach for each crime type fails to recognize this correlation and may ignore information shared in those missing variables.In ecological analysis, ignoring correlation structures between outcomes may distort the estimates of regression coefficients and result in biased and inefficient parameter inferences [19].Nonetheless, to the best of our knowledge, there has been little work on multivariate Bayesian crime modeling in criminological research.
This paper uses a multivariate Bayesian spatial modeling approach to jointly model the counts of two types of crime, i.e., burglary and non-motor vehicle theft, and explore the geographic pattern of crime risks at a small-area scale.Burglaries and non-motor vehicle thefts can be interrelated as they belong to the same theft category in China and may share some omitted variables.The multivariate Bayesian spatial model takes into consideration the correlation between these two types of crime.This is of both theoretical and practical importance.From the theoretical perspective, there are gains in precision and efficiency of estimates through this joint modeling method achieved by pooling all the available data from different crime sources.Multivariate Bayesian models can take advantage of correlations across different types of crime.In addition, the multivariate spatial model provides information not only on similarities but also on differences in the effects of covariates on crimes.The model could also be used to identify crime-specific risk factors.From a practical perspective, both burglaries and non-motor vehicle thefts are common and have a high incidence in China.Police departments and the public are greatly concerned about these two types of crime and are interested in identifying relevant risk factors.The information gained from the joint analyses could help develop a greater understanding of the spatial distribution and determinants of crimes, contributing to more efficient and effective crime intervention and prevention strategies.Through the joint multivariate crime modeling approach, these issues can be addressed in a unified framework.
The remainder of this manuscript is organized in the following manner.First, we briefly describe the study region and relevant data.Next, we explain the multivariate Bayesian spatial modeling approach in detail.Third, results of the analysis are shown, followed by a discussion of the results.The strengths and limitations of the study, as well as directions for feature research are concluded in the final section.

Study Region
This study was carried out in Jianghan District, Wuhan, China.Wuhan, the largest city in Central China, is an important industrial base and transportation hub.The Jianghan District of Wuhan lies at the intersection of the Yangtze River and the Han River.It is the most prosperous and populated urban district in the city.The neighborhood was the unit of analysis in the current research.It is a direct sublevel of a sub-district, which is one of the smallest political divisions in China.The Jianghan District consists of 116 neighborhoods; two neighborhoods were excluded from the study because of missing data.

Data
Data on burglaries and non-motor vehicle thefts came from the city's 110-reporting system managed by Wuhan Public Security Bureau.The 110-reporting system is the primary source of official crime information in China, as 110 is the emergency telephone number for the public to report crimes and emergencies to police.A total of 1747 burglaries and 2341 non-motor vehicle thefts occurred in the study region from 1 January 2015 to 31 December 2015, respectively.Both types of incidents were unevenly distributed among neighborhoods (see Table 1).We used these incident data as dependent variables.Independent variables were selected based on social disorganization theory [20] and routine activity theory [21].Social disorganization theory suggests that crime levels in neighborhoods are strongly correlated with local ecological characteristics.Socioeconomic stress can decrease social organization.Socially disorganized neighborhoods have lower levels of social control, which will decrease the monitoring of criminal activity and further lead to occurrences of crime [22,23].Routine activity theory identifies three elements necessary for a crime: motivated offenders, suitable targets, and the absence of capable guardians [24,25].If any of these is missing, a crime will be less likely to occur.According to these two theories, six variables at the neighborhood level were assumed to be potential risk factors.They were population density, unemployment rate, percentage of highly educated, percentage of young males, policing and bar density.
Population density was measured by the number of residents per square kilometer.It was included to control for differences in neighborhood sizes and populations [26,27].It is generally believed that a larger population density usually leads to a higher crime rate [28,29].However, higher population density also increases street monitoring by residents [30], deterring anyone from committing crime.
Unemployment rate was calculated as the number of unemployed residents between ages 18-60 divided by the number of all the residents within the same age range in the neighborhood.It is usually believed that unemployment rate is positively related to crime and researches have used it as one indicator of socioeconomic status [14,31,32].Percentage of highly educated was the percentage of population age 20 or older with a high school diploma or equivalent.This variable is also often used to measure economic deprivation [14,32,33].Therefore, in our analysis, the unemployment rate and percentage of highly educated were used as a proxy for socioeconomic status.Neighborhoods with a high unemployment rate and a low percentage of highly educated inhabitants, relatively speaking, likely have a lower socioeconomic status.According to social disorganization theory [20], there is often a positive correlation between low socioeconomic status and crime.
Percentage of young males referred to the percentage of males between ages 15-24.This factor was considered because previous research has shown that young males are hugely overrepresented among those engaged in criminal activities [23,34,35].
Policing was measured by the number of community policing and monitoring rooms in each neighborhood per 10,000 people.Policing rooms are set up by local police stations and staffed by several police officers to guarantee neighborhood's security.The relationship between policing and crime has been investigated in previous studies [15,27].
Bar density was defined as the number of bars per 10,000 people in each neighborhood.We included this land-use variable in the analysis to test if it could explain the spatial pattern of crime risks.Researchers have long been interested in land-use variables in crime analysis [36,37].
Descriptive statistics, including minimum, maximum, mean and standard deviation (SD) of the dependent and independent variables at the neighborhood level are shown in Table 1.
Table 2 summarizes the correlation matrix of the six independent variables.This table indicates that the pair of population density and unemployment rate had the most significant correlation coefficient (0.420).According to Evans [38], this correlation level is moderate while those of other pairs of variables are weak or very weak.Therefore, we did not consider the multicollinearity issue and used all six variables in the modeling.

Modeling Framework
The dependent variable was the number of crimes by type, assumed to follow a Poisson distribution: where k refers to the type of crime, Y ik is the observed number of crimes of the type k in neighborhood i, E ik r ik is the expected value of the Poisson distribution, E ik is the number of crimes that could be expected, and r ik is the corresponding area-specific and type-specific crime risk.In our Bayesian spatial modeling approach, the logarithm of the expected value was linked to the independent variables: where log(E ik ) is an offset term for crime type k in neighborhood i with the regression coefficient being one, β k0 is the intercept for crime type k, X i1 , • • • , X il are the observations of the independent variables in neighborhood i, β k1 , • • • , β kl are the corresponding regression coefficients for crime type k, U ik is an unstructured random effects term to account for the overdispersion that may arise when small-area count data are modeled, and S ik is a spatially structured random effects term to account for spatial autocorrelation.

Prior Specification for the Univariate Model
In order to fully express the multivariate Bayesian spatial model, we describe prior specifications used in the univariate model and contrast them to their multivariate counterparts.Prior specification for unknown parameters is important in Bayesian inference as posterior distributions of parameters are deduced by combing prior knowledge and data.β k0 was assigned a uniform prior dflat() due to a sum-to-zero constraint on the random effects [39]: β kl , the regression coefficient, was specified to follow a non-informative normal distribution with a mean of zero and a variance of 10,000 as we had no genuine prior expectations for direction and magnitude of the independent variable's effect: The univariate Bayesian spatial model ignores correlations across different types of crime as it assumes that the random effects for different types of crime are independent.In the case of the univariate model, U ik was specified as an independent normal distribution with an expected mean of zero and a variance of σ 2 uk : S ik , the spatially structured random effects term, was determined through an intrinsic conditional autoregressive (ICAR) distribution that is widely used in spatial statistics [40,41].Under the ICAR specification, this will result in: where S ik = ∑ j =i ω i,j S jk /n i , ω i,j = 1 if i and j are adjacent and otherwise 0, and n i is the number of neighborhoods adjacent to neighborhood i.Specifically, according to the specification of the ICAR model, the conditional distribution of S ik given the remaining components (S −ik ) is normal with mean S ik and variance σ 2 sk /n i .The variation of S ik is controlled by the overall variance parameter σ 2 sk .Regarding hyperparameters σ 2 uk and σ 2 sk , the reciprocals of these two parameters were assumed to follow a commonly used Gamma distribution [42]: The WinBUGS code for the univariate model can be found in the Supplementary Materials.

Prior Specification for the Multivariate Model
The specifications from Equation (4) to Equation (9) are essentially equivalent to the specifications of the univariate Poisson spatial model for each crime type, i.e., burglary and non-motor vehicle theft, respectively.However, these specifications provide the possibility of being directly comparable with the multivariate joint specifications.
The prior specifications of β k0 and β kl for the multivariate Bayesian spatial model were the same as those for the univariate model.In contrast to the univariate model assuming that the random effects for different types of crime are independent, the multivariate model takes account of correlations across random effects.Therefore, the main differences between the multivariate model and the univariate model are the prior specifications of the random effects.For the multivariate spatial model, U i = (U i1 , U i2 ) T was assumed to follow a multivariate normal distribution: (10) where µ i is a vector with all elements zeroes µ i = (0, 0) T and ∑ u is the covariance matrix estimated by a Wishart distribution: where Ω = 0.01 0 0 0.01 and d = 2 are the scale matrix and the degrees of freedom, respectively.
For computational convenience, the Wishart distribution is commonly used as it is the conjugate prior for the inverse of the covariance parameters of multivariate normal distributions [43,44].The scale matrix, Ω, is often set to be a scaled identify matrix with a scaling factor (see [45][46][47], for example).
In our analysis, we chose 0.01 as the scaling factor, i.e., setting Ω 's diagonal entries to Ω ii = 0.01, a specification adopted to run sensitivity analysis for the precision matrix of a multivariate normal distribution [45].The parameter d, the degrees of freedom, was set equal to 2 in order to make the prior on ∑ u −1 minimally informative [48,49].
As to the spatially structured random effects S i = (S i1 , S i2 ) T , the ICAR specification in the univariate model extends naturally to a multivariate ICAR specification by replacing the univariate normal conditional distribution with a multivariate conditional distribution: where S i is the mean vector S i = (∑ j =i ω i,j S j1 /n i , ∑ j =i ω i,j S j2 /n i ) T , and ∑ S is the 2 × 2 covariance matrix, respectively.As in the univariate case, ω i,j = 1 if i and j are adjacent and 0 otherwise, and n i is the number of neighbors of neighborhood i. Diagonal elements of ∑ S , i.e., ∑ 11 and ∑ 22 , are the conditional variances of S i1 and S i2 respectively, and off-diagonal elements, i.e., ∑ 12 and ∑ 21 , represent the conditional within-neighborhood covariance between S i1 and S i2 .This specification of S i gives the multivariate equivalent of the standard ICAR specification [40].The multivariate ICAR model is building on multivariate conditional autoregressive (MCAR) models described by Mardia [50].It is the intrinsic version of the Multivariate Conditional Autoregressive (MCAR) model [51] and has been used in spatial analysis [46,52].Similar to the specification of ∑ u , we adopted the same approach to determine ∑ S ; a Wishart distribution was chosen for its inverse ∑ s −1 [45,48,49]: where R = 0.01 0 0 0.01 and d = 2 are the scale matrix and the degrees of freedom, respectively.
The WinBUGS code for the multivariate model can be found in the Supplementary Materials.

Model Implementation and Assessment
Models were fitted by the MCMC simulation approach in WinBUGS 1.4.3.Before they were fitted, all independent variables were standardized (centered around the mean and then divided by the standard deviation) to improve convergence [51].For each model, two parallel sampling chains with different initial values were run.We evaluated the convergence of each model by visual inspection of the trace plots and the Gelman-Rubin convergence statistic [53].Convergence occurred after the first 50,000 iterations.Each chain was run for a further 100,000 iterations, generating 200,000 samples with acceptable Monte Carlo errors (<5% of the posterior standard deviation) for posterior estimations.We assessed different models using the deviance information criterion (DIC) [54].Broadly, DIC is a generalization of Akaike Information Criterion [55] and it is defined as: where D is the posterior mean of the deviance and p d is the number of effective parameters in the model.DIC takes both model fit (summarized by D) and model complexity (captured by p d ) into consideration when comparing models.Smaller values of DIC suggest better-fitting models.

Comparison of the DICs
Table 3 presents the results of the univariate and the multivariate regression models fitted in WinBUGS.Both models included an unstructured random effects term and a spatially structured random effects term.The major difference is the way in which the random effects were specified; one is in the univariate form, while the other is in the multivariate form.The univariate model had a DIC of 1344.709; when extended to the multivariate model, the DIC dropped to 1325.992, a difference of 18.717.This difference in DICs can be explained by the posterior deviance D and the number of effective parameters p d .The D for the multivariate model was 1145.07, as compared to 1150.30 for the univariate model.The higher D for the univariate model indicates a poor fit relative to the fit of the multivariate model.The p d in the univariate model was 194.409, while extending the univariate random effects to the multivariate form made p d 13.487 points lower.This reduction in p d indicates that the random effects might be correlated.Indeed, once the correlation between random effects was considered in the multivariate model, the p d dropped, largely explaining the improvement in DIC.In short, the multivariate model was more parsimonious than the univariate alternative.A model with a better fit and fewer parameters is always preferred.Therefore, in terms of DIC, the overall model comparison criterion that takes into account both model fit and model complexity, the multivariate model was better suited to our data.We studied the sensitivity of the results to the choice of hyperparameters.Regarding the univariate model, we note that standard deviation uniform priors are now commonly used [56].Thus, we used a non-informative uniform hyperprior distribution dunif(0,100) for the standard deviation parameters of U ik and S ik [16].It yielded a DIC with a difference of 5.077 (1349.786versus 1344.709).For the multivariate model, we selected the two-dimensional identity matrix, an uninformative prior often used in practice, for Ω and R [45,47].It yielded a DIC with a difference of 4.198 (1330.190 versus 1325.992).The sensitivity analysis suggests that the results of the univariate and the multivariate models were both not overly sensitive to the selected hyperprior specifications.

Analysis of Risk Factors
Regarding the independent variables, population density and bar density had a clear association with both types of crime at the 95% credible interval (CI) in both the univariate model and the multivariate model.The other four variables, i.e., the unemployment rate, percentage of highly educated population, percentage of young males and policing, were not significant for the two types of crime in both models.The univariate model and the multivariate model identified the same set of significant variables for both crimes.However, the precisions of all the regression coefficients from the univariate model were overestimated as compared to those from the multivariate model.Overall, the regression coefficients were much changed after extending the univariate random effects to their multivariate counterparts.This indicates that the correlation between burglaries and non-motor vehicle thefts had a great effect on coefficient estimates.Since the multivariate model better fit the data, our analyses are based on this model.
Population density had a negative effect on both types of crime.Bar density was found to be positively correlated with both burglaries and non-motor vehicle thefts at the 95% CI.Generally, holding all other variables constant, the relative risk of burglary and non-motor vehicle theft will decrease by 28.12% and 41.32%, respectively, with an increase of one standard deviation in population density (Equations ( 15) and ( 16)).exp(−0.3302)− 1 = −0.2812(15) exp(−0.5331)− 1 = −0.4132(16) Using the same calculation method, we estimated that the relative risk of burglary and non-motor vehicle theft will increase by 67.30% and 46.96%, with an increase of one standard deviation in bar density, respectively (Equations ( 17) and ( 18)).exp(0.5146)− 1 = 0.6730 (17) exp(0.3850)− 1 = 0.4696 (18)

Correlation between the Random Effects
Table 4 shows the posterior means of the covariance/correlation matrix of the random effects.Although the correlation between the unstructured random effects was low (0.1398), the correlation between the spatially structured random effects was high (0.7664).Possibly due to the dominance of the spatially structured random effects, the posterior correlation between the total random effects was also high (0.7953), indicating the close relationship between the two types of crime in the data.As a consequence, like any statistical model, this correlation needs to be incorporated into the estimation.In the multivariate model, by calculating the inner product of the two crime risks, the mean and the standard deviation of type-specific crime risks, we further computed the correlation between the crime risks.The posterior mean of the correlation was 0.7592, again suggesting a strong shared geographical pattern of risks between burglaries and non-motor vehicle thefts.

Mapping the Relative Risks
The posterior mean risks of burglary and non-motor vehicle theft estimated by the multivariate model are shown in Figures 1a and 1b, respectively.In contrast to the standardized crime ratio (SCR), a measure of relative risks obtained by dividing the observed number of crime incidents to the expected number, the risks shown in Figure 1 were estimated while accounting for the effects of the independent variables, the unstructured random effects, and the spatially structured random effects.
No neighborhood was assigned a risk value of exactly zero in Figure 1, although the SCRs of some neighborhoods would be zero as no crimes occurred in these neighborhoods.This difference between the SCRs and the risks shown in Figure 1 was mainly due to the spatial smoothing effect from the spatially structured multivariate random effects.The multivariate model can borrow strength not only from the estimates of the other crime type in the same neighborhood but also through the spatial correlation, from those in the adjacent neighborhoods.This "borrowing strength" effect helps stabilize risk estimations.In addition to the posterior means, we also mapped the 2.5% percentile of the posterior distribution of s (Figure 2).Percentiles of the posterior distribution of parameters of interest are useful in interpreting the results of a Bayesian analysis.In the current study, the 2.5% percentile provided valuable information about the pattern of crime risks.The different shadings shown in Figure 2 refer to the value of such that 97.5% of its simulated values were greater than the specified value.Figure 2a shows that there were 18 neighborhoods with a posterior 2.5% percentile value of burglary risk great than 1.0.Of particular note is that the percentile value in two neighborhoods was higher than 10.0. Figure 2b indicates that there were 23 neighborhoods having the 2.5% percentile value of non-motor vehicle theft risk greater than 1.0, and one neighborhood having a percentile value greater than 10.0.In short, nine neighborhoods had the 2.5% percentile values of both burglary risk and non-motor vehicle theft risk greater than 1.0 and one neighborhood had percentile values higher than 10.0 for both types of crime.These neighborhoods need to be further investigated by researchers and should be given adequate attention by the relevant authorities.In terms of the type-specific relative risks, Figure 1a indicates that three neighborhoods had a burglary risk value higher than 10.0 and 38 neighborhoods had a risk value greater than 1.0.With regard to non-motor vehicle theft, there were also three neighborhoods with a risk value higher than 10.0 as can be seen from Figure 1b.In addition, 43 neighborhoods had a non-motor vehicle theft risk greater than 1.0.Taken together, two neighborhoods had both crime risks higher than 10.0.About 25.44% of all the neighborhoods (29 out of 114 neighborhoods) had both crime risks greater than 1.0, indicating they had crime risks higher than the average for the whole study region.
In addition to the posterior means, we also mapped the 2.5% percentile of the posterior distribution of r ik s (Figure 2).Percentiles of the posterior distribution of parameters of interest are useful in interpreting the results of a Bayesian analysis.In the current study, the 2.5% percentile provided valuable information about the pattern of crime risks.The different shadings shown in Figure 2 refer to the value of r ik such that 97.5% of its simulated values were greater than the specified value.Figure 2a shows that there were 18 neighborhoods with a posterior 2.5% percentile value of burglary risk great than 1.0.Of particular note is that the percentile value in two neighborhoods was higher than 10.0. Figure 2b indicates that there were 23 neighborhoods having the 2.5% percentile value of non-motor vehicle theft risk greater than 1.0, and one neighborhood having a percentile value greater than 10.0.In short, nine neighborhoods had the 2.5% percentile values of both burglary risk and non-motor vehicle theft risk greater than 1.0 and one neighborhood had percentile values higher than 10.0 for both types of crime.These neighborhoods need to be further investigated by researchers and should be given adequate attention by the relevant authorities.
neighborhoods was higher than 10.0. Figure 2b indicates that there were 23 neighborhoods having the 2.5% percentile value of non-motor vehicle theft risk greater than 1.0, and one neighborhood having a percentile value greater than 10.0.In short, nine neighborhoods had the 2.5% percentile values of both burglary risk and non-motor vehicle theft risk greater than 1.0 and one neighborhood had percentile values higher than 10.0 for both types of crime.These neighborhoods need to be further investigated by researchers and should be given adequate attention by the relevant authorities.

Discussion
This research applied a multivariate Bayesian spatial modeling approach to jointly model the counts of burglaries and non-motor vehicle thefts and estimate the relative risks of these two types of crime at a small-area scale.The multivariate model accounts for the correlation between different types of crime.This correlation cannot be recognized when separate univariate models are estimated for each crime type.Although the univariate model identified the same set of significant covariates for both crimes as the multivariate model did, a smaller DIC for the multivariate spatial model suggests that it was better supported by our data (Table 3).In addition, the precisions of all the regression coefficients from the univariate model were overestimated as compared to those from the multivariate model (Table 3).Therefore, the multivariate model is superior to the univariate model.
Our analyses show that population density and bar density had a clear association with both types of crime in both the univariate and the multivariate models (Table 3).Population density was negatively correlated with both crime risks.The clear association between population density and crime found in the present study is in accordance with other studies showing that residential density can lower crime rates [26,27].This may be due to the closer and easier street monitoring possible when there is a larger population density [30].Burglary and non-motor vehicle theft risks were higher in neighborhoods that had more than the average bar density.The linkage between bars and crimes is generally supported in the literature [57].The unemployment rate and percentage of highly educated were not significantly related to both types of crime in either model.This is somewhat unexpected as it is not in accordance with social disorganization theory.One possible reason may be that these two variables cannot fully represent socioeconomic status.Issues regarding the ecological fallacy [58,59] and the modifiable area unit problem [60] also warrant consideration.However, similar results have been found in an analysis of other crime data [32].Policing did not have a clear association with both types of crime.This is not in line with Faria et al. [27], who found that there was a positive relation between policing and burglaries.Nevertheless, our results are not entirely unexpected.Policing can prevent crime but a larger policing presence is typically a reaction to high crime areas.Therefore, there may not be a clear association between policing and crime.
Through the joint modeling of burglaries and non-motor vehicle thefts, our results suggest that there was a tight relationship between the two types of crime (Table 4).The correlation between the relative risks of the two crimes was high (0.7592), indicating a strong shared geographical pattern of risks between burglaries and non-motor vehicle thefts.
We mapped the means (Figure 1) and the 2.5% percentile (Figure 2) of the posterior distribution of crime risks estimated by the multivariate Bayesian spatial model.By demonstrating the spatial distribution of crime risks and highlighting areas with high crime risks, the maps shown in Figures 1 and 2 provide a scientific basis for policy formulation and resource allocation.
We note that some other models, such as the shared component model and the smoothed analysis of variance (SANOVA) model, may also be used to jointly model different types of crime.The shared component model can separate the underlying risk surface for each crime into a shared and a crime-specific component [61].However, when the shared component model is extended to the joint analysis of more than two types of crime, the number of possible permutations of shared and specific components may rapidly become prohibitive [61].Moreover, the interpretation of a shared component in such cases is more difficult [51].The SANOVA model is also capable of modeling a multivariate outcome [62][63][64].According to this model, a series of orthogonal contrasts (in the sense of analysis of variance) is considered [63].The choice of these contrasts is a main problem with the SANOVA model, as the degree of benefit obtainable through fitting the model mainly depends upon how appropriate this choice is, and often there are no clear standards and criteria to help choose them [64].The multivariate ICAR model adopted in our analysis is a natural extension of the commonly used ICAR model in univariate studies.It can easily be extended to the joint analysis of three or more types of crime.

Conclusions
Our study represents a contribution to the existing criminological literature as the multivariate Bayesian spatial modeling approach is seldom used in crime analysis.In contrast to the univariate model, the multivariate Bayesian spatial model not only takes into account correlations across different types of crime in the same spatial unit but also takes advantage of its neighbors to borrow strength.This "borrowing strength" effect can overcome the small number problem and contributes to stable risk estimations.
The results of our study have practical implications.The use of the multivariate Bayesian spatial modeling approach provides a new perspective to better understand the risk factors associated with different types of crime (Table 3) as well as the relation between crimes (Table 4).It can highlight areas with high crime risks (Figures 1 and 2) and therefore provides a scientific basis for crime prevention and intervention strategies.Generally, police in charge of burglaries and non-motor vehicle thefts in the Jianghan District should focus on neighborhoods that have a high bar density.In addition, neighborhoods showing the highest crime risks should be given adequate attention more generally.A further investigation, identifying environmental exposures specific to these parts of the Jianghan District showing a high crime risk, may be worthwhile.This helps identify other features of the neighborhood context to explain high levels of risk and thus can be used for a better targeted prevention and intervention efforts.Additionally, as our analysis was conducted at a small-area scale, it can inform an efficient use of resources.Rather than allocate resources through larger areas such as sub-districts, understanding the distribution of crime risks at the neighborhood level will facilitate more precise targeting of resources such that they are used most efficiently and effectively.
This study also has some limitations.First, our crime data were obtained from the city's 110-reporting system.Flaws in the data, such as data entry errors and the issue of missing data related to non-reporting of crimes, have been acknowledged [65].However, despite these limitations, the 110-reporting system does provide the most authoritative source of crime information in China.Second, we used residential population as the at-risk population to calculate expected crime counts.We did not take into account daytime changes in population size due to commuting.This may influence the risk estimations and bias the assessment of contextual influences on crime risks in ecological crime studies as indicated by recent research [66][67][68][69].Third, the measures used to represent socioeconomic status (i.e., the unemployment rate and percentage of highly educated) may not be powerful enough, as other related variables, such as average household income [27], percentage of individuals below the poverty line [70] and percentage of high occupational status [14], were not included.In addition, socioeconomic status is only one indicator of social disorganization.Other indicators commonly used to operationalize social disorganization include ethnic heterogeneity, family disruption, socioeconomic status and residential mobility [4,33].However, due to unavailability of data, these indicators are not included in our analysis.Another limitation of the study is that we did not account for the role of the spatial configuration in explaining crime risks.Recent research has shown that spatial configuration can be a substantial factor in creating crime risks and in revealing the sources of spatial autocorrelation in crime data [71].
Despite the limitations, we would like to emphasize that the main purpose of our paper is to present and illustrate a multivariate Bayesian spatial model.Previous crime data are often analyzed by modeling each crime type separately without accounting for correlations that may exist among different types of crime.The multivariate model presented in the current research can cope with both the correlation structure and overdispersion in the crime data.It acknowledges dependence between different types of crime as well as dependence across space.In addition, this model can easily be extended to jointly analyze more than two types of crime.
There are several important avenues for future work.First, additional tests of the multivariate model with different crime data are recommended.Second, the current multivariate model should be extended to incorporate more than two types of crime to deepen our understanding of crime patterns and the correlations among different crimes.Third, future research should fully investigate if spatial configuration is a relevant variable in spatial models of urban crime as this is an often under evaluated topic in crime data modeling.Fourth, as shown by Helbich and Jokar Arsanjani [72], crime patterns change over time.It is advisable to collect more years of crime data and extend the multivariate spatial model into a spatio-temporal one to analyze the change of crime patterns over years.

15 Figure 1 .
Figure 1.The posterior mean risks of (a) burglary and (b) non-motor vehicle theft estimated by the multivariate model.

Figure 1 .
Figure 1.The posterior mean risks of (a) burglary and (b) non-motor vehicle theft estimated by the multivariate model.

Figure 2 .
Figure 2. The 2.5% percentile of the posterior distribution of (a) burglary risks and (b) non-motor vehicle theft risks estimated by the multivariate model.

Figure 2 .
Figure 2. The 2.5% percentile of the posterior distribution of (a) burglary risks and (b) non-motor vehicle theft risks estimated by the multivariate model.

Table 1 .
Descriptive statistics of the dependent and independent variables.

Table 2 .
Correlation matrix of the independent variables.

Table 3 .
Results of the univariate and the multivariate models.

Table 4 .
Posterior covariance/correlation matrix of the random effects: upper triangle, covariance of the random effects; lower triangle, correlation of the random effects.