Spatiotemporal Suicide Risk in Germany: A Longitudinal Study 2007–11

Despite comprehensive prevention programs in Germany, suicide has been on the rise again since 2007. The underlying reasons and spatiotemporal risk patterns are poorly understood. We assessed the spatiotemporal risk of suicide per district attributable to multiple risk and protective factors longitudinally for the period 2007–11. Bayesian space–time regression models were fitted. The nationwide temporal trend showed an increase in relative risk (RR) of dying from suicide (RR 1.008, 95% credibility intervals (CI) 1.001–1.016), whereas district-specific deviations from the grand trend occurred. Striking patterns of amplified risk emerged in southern Germany. While the number of general practitioners was positively related (RR 1.003, 95% CI 1.000–1.006), income was negatively and non-linearly related with suicide risk, as was population density. Unemployment was associated and showed a marked nonlinearity. Neither depression prevalence nor mental health service supply were related. The findings are vital for the implementation of future suicide prevention programs. Concentrating preventive efforts on vulnerable areas of excess risk is recommended.


Specification of space-time regressions
To model contributing risk and protective factors from 2007 to 2011 and to investigate spatiotemporal suicide risk, hierarchical Bayesian models were implemented 1,2 . The application was restricted to the parametric time trend model 3 and the non-parametric dynamic counterpart 1,4 which were found to be superior to more complex models 5,6 .
As suicide is a count, it is valid to assume that these counts arise from a Poisson distribution. Let y it be the observed suicide cases in area i (i=1,…,402) at time t (t=2007,…,2011), it ρ denotes a rate, and it E represent the expected number of cases 1,2 . Then the implemented mode is expressed as (equation Model 1a follows Bernardinelli et al. (1995) who proposed the following spatiotemporal model with a parametric linear time trend (equation 2): where the intercept α represents the area-wide relative risk and k β refers to the regression coefficient of covariate k. Whereas the covariate income, unemployment rate, and population density is timevarying, the remaining covariates are kept temporally constant. In model 1b, significantly associated linear covariates are replaced through second-order random walks to model non-linear effects 7, 8 . To obviate violation of the model assumption of spatial independence, it was of paramount importance to explicitly model area-specific spatial effects; that is, spatially adjacent districts have an associated risk.
Here, following Besag et al. (1991), i v is a spatially structured residual effect for each district modeled as intrinsic conditional autoregressive specification. An unstructured residual effect i u 2 models spatially uncorrelated extra variability provoked by unobserved aspatial variables and is assumed to follow a Gaussian distribution 10 . The smoothed district-specific relative risk compared to entire Germany is ) exp( Districts are neighbours when a common boundary is shared. Rather than pooling data over time, models 1a and 1b consider temporal dependencies. ψ represents the grand (i.e., study area-wide) temporal trend, whereas the differential component, i δ , captures district-related temporal deviations modelled as independent and Gaussian distributed random variable. Interpretatively, a differential component below 0 indicates a less pronounced trend than the grand trend, and a value of above 0 refers to a more pronounced trend 2 .
Model 2 relaxes the linearity assumption of the temporal effect in models 1a and 1b by means of a nonparametric dynamic time trend 4 (equation 3): Here, the same parametrization applies as above except that t γ now refers to a temporally structured effect modelled as second-order random walk and t φ to a temporally unstructured effect. The temporally structured effect aims to resemble the trend of adjacent districts.
Bayesian inference was carried out with the integrated nested Laplace approximation, which is a highly accurate and computationally fast alternative to Markov chain Monte Carlo methods 11 .
Although the default prior distributions were specified for the model parameters, the hyperparameters for the spatial effects and random walks were scaled to achieve a less ad hoc selection 12 . For all models, relative risk estimates were obtained based on exponentiating their posterior means together with the 95% credibility intervals (CI) to obtain parameter uncertainty. In order to assess the quality of the models, three statistics were utilized: a) The deviance information criterion (DIC) to judge the goodness-of-fit and to compare candidate models 13