Analysis of crime data in Austria

I looked at police-recorded crime in Austria using Eurostat data: regional counts for 2021, the mix of offence types across regions, and a few hypothesis tests.

data preprocessing

I used the NUTS 3 regions to split Austria into smaller areas.

# Load necessary libraries
library(eurostat)
library(ggplot2)
library(psych)

id = 'crim_gen_reg'
crim_data = get_eurostat(id=id)

# Filter for data from the year 2021
data_2021 = subset(crim_data, format(TIME_PERIOD, '%Y') == '2021')

# Filter for Austrian NUTS 3 regions
at_data = data_2021[grepl('^AT[0-9]{3}$', data_2021$geo), ]
df = subset(at_data, select = c(unit, iccs, geo, values))
df = label_eurostat(df)

# The subcategories 'Burglary of private residential premises' and 
# 'Theft of a motorized land vehicle' are already included in 'Burglary' and 'Theft'.
# We will exclude them to avoid duplication.
df = subset(df, !(iccs %in% c('Burglary of private residential premises',
                              'Theft of a motorized land vehicle')))

# Separate data into absolute numbers and per 100k inhabitants
nr_df = subset(df, df$unit == 'Number', select = c(iccs, geo, values))
pht_df = subset(df, df$unit == 'Per hundred thousand inhabitants', select = c(iccs, geo, values))

# Aggregate data by crime category (iccs) and region (geo)
nr_iccs_df = aggregate(list(values = nr_df$values), list(iccs = nr_df$iccs), sum)
nr_geo_df = aggregate(list(values = nr_df$values), list(geo = nr_df$geo), sum)

pht_iccs_df = aggregate(list(values = pht_df$values), list(iccs = pht_df$iccs), mean)
pht_geo_df = aggregate(list(values = pht_df$values), list(geo = pht_df$geo), mean)

initial data exploration

First, the 2021 counts by offence type and region.

At this level, the regions are groups of political districts (Gruppen von Politischen Bezirken), as shown on the map above.

categories of criminal offenses

From 2008 onwards, the statistics include police-recorded offences for homicide, assault, sexual violence, robbery, burglary, (of which) burglary of residential premises, theft, (of which) theft of motorized land vehicle. [src]

rows = nr_iccs_df[order(-nr_iccs_df$values),]
row.names(rows) = 1:5
rows

theme_set(theme_gray(base_size = 14))
options(repr.plot.width=14, repr.plot.height=6)
ggplot(nr_iccs_df, aes(x=reorder(iccs, values), y=values)) + ggtitle('Categories of criminal offences') +
    geom_bar(stat='identity') + ylab('count') +
    geom_text(aes(label=values), hjust=-0.3) +
    coord_flip() + theme(axis.title.y = element_blank()) +
    expand_limits(y = c(0, 75000))
ggplot(pht_df, aes(x=reorder(iccs, values), y=values, fill=iccs)) +
    geom_boxplot(outlier.color='red', show.legend=F) + ylab('relative count [per 10^5 inh.]') +
    coord_flip() + theme(axis.title.y = element_blank()) +
    geom_jitter(color='black', size=0.1, alpha=0.8, show.legend=F)

Theft is by far the most common offence, with 73,213 cases against 59 intentional homicides. To compare regions of different size, I also plotted the rates per 100,000 inhabitants. Theft rates still vary considerably between regions, with Wien and Linz-Wels as outliers.

regions by NUTS 3 division

rows = merge(x = nr_geo_df, y = pht_geo_df, by = 'geo')
rows = rows[order(-rows$values.x),]
colnames(rows) <- c('geo', 'count', 'relative.count')
row.names(rows) = 1:35

cat('Top 5')
head(rows, 5)

cat('Bottom 5')
tail(rows, 5)

summary(nr_geo_df)
r = describe(nr_geo_df$values, skew=F, IQR=T, ranges=F)
rownames(r) = c('values')
r$var = c(var(nr_geo_df$values))
r[,c(2,4,7,6)]

ggplot(rows, aes(x=count/relative.count, y=count)) + geom_point() + ggtitle('Inhabitants v. Crime') +
    xlab('inhabitants [10^5]') + ylab('crime count')
ggplot(rows[rows$geo != 'Wien', ], aes(x=count/relative.count, y=count)) + geom_point() + ggtitle('Inhabitants v. Crime (w/o Wien)') +
    xlab('inhabitants [10^5]') + ylab('crime count')

Top 5

Bottom 5

          geo            values       
 Length:35          Min.   :  112  
 Class :character   1st Qu.: 1034  
 Mode  :character   Median : 1970  
                    Mean   : 4287  
                    3rd Qu.: 3352  
                    Max.   :63198  

The largest counts are in Wien (63,198), Linz-Wels (11,245), and Graz (8,469). At the other end are Mittelburgenland (164) and Lungau (112). Much of this difference is likely explained by population, which I test in hypothesis 1 below.

rows = merge(x = nr_geo_df, y = pht_geo_df, by = 'geo')
rows = rows[order(-rows$values.y),]
colnames(rows) <- c('geo', 'count', 'relative.count')
row.names(rows) = 1:35
head(rows, 5)

options(repr.plot.height=6)
ggplot(pht_geo_df, aes(x=values)) +
    geom_histogram(bins=8) + xlab('Relative crime count') + ylab('Frequency') +
    geom_rug(aes(values, y = NULL), length = unit(0.02, "npc")) +
    geom_boxplot(outlier.color='red', show.legend=F, position = position_nudge(y = -0.2))

Per capita, the ranking changes. Bludenz-Bregenzer Wald and Salzburg und Umgebung have relatively low total counts but rank second and third after Vienna.

analyzing the relationship between region and crime type

ct = xtabs(formula=values ~ geo + iccs, data=nr_df)
addmargins(ct)

The contingency table mostly reflects region size and how common each offence is. The relative mix of offence types looks similar across regions, so I tested whether it actually is.

Fisher’s exact test

I used Fisher’s exact test at the 5% level to test whether the offence type is independent of the region. I chose it over the χ² test because some cells have zero counts.

\(H_0\): Each row (region) is a realization of the same distribution (crime category), i.e., \(p_{ij} = p_{i:} \cdot p_{:j}\) for \(i=1,\ldots,35\) and \(j=1,\ldots, 5\).

\(H_A\): \(H_0\) is not true.

fish = fisher.test(ct, simulate.p.value = TRUE)
fish

  Fisher's Exact Test for Count Data with simulated p-value (based on
  2000 replicates)

data:  ct
p-value = 0.0004998
alternative hypothesis: two.sided

The test rejects the null hypothesis at 5%: the distribution of offence types differs between regions, even though this is not obvious from the table. With this many observations, even small differences are significant, so the test says little about how large the differences are.

hypothesis testing

hypothesis 1: correlation between population and crime count

ggplot(rows, aes(x=count/relative.count, y=count)) + geom_point() + ggtitle('Inhabitants v. Crime') +
    xlab('inhabitants [10^5]') + ylab('crime count')
ggplot(rows[rows$geo != 'Wien', ], aes(x=count/relative.count, y=count)) + geom_point() + ggtitle('Inhabitants v. Crime (w/o Wien)') +
    xlab('inhabitants [10^5]') + ylab('crime count')

The plots above suggest that more populous regions have more recorded crime. I tested this with Spearman’s rank correlation at the 5% level.

\(H_0 :\) There is zero correlation between the number of inhabitants and the frequency of crime in a region: \(\rho_S = 0\)

\(H_A :\) There is a positive correlation between the number of inhabitants and the frequency of crime in a region: \(\rho_S \> 0\)

x = rows$count/rows$relative.count
y = rows$count

round(cor(x, y, method='spearman'), 4)
cor.test(x, y, method='spearman', alternative='greater')

0.8641

  Spearman's rank correlation rho

data:  x and y
S = 970, p-value = 1.586e-08
alternative hypothesis: true rho is greater than 0
sample estimates:
       rho 
0.8641457 

The Spearman correlation is 0.86, a strong positive relationship. The test rejects zero correlation at 5%, consistent with the plots.

hypothesis 2: comparing homicide counts in Austria and Slovakia

cd = crim_data[grepl('^AT[0-9]{3}$', crim_data$geo), ]
cd = label_eurostat(cd, fix_duplicated = TRUE)
cd = subset(cd, unit == 'Number')
cd = subset(cd, freq == 'Annual')
cd = subset(cd, iccs == 'Intentional homicide')
cd = aggregate(list(values = cd$values), list(TIME_PERIOD = cd$TIME_PERIOD), sum)
mean.AT = mean(cd$values)
median.AT = median(cd$values)
colnames(cd) = c('year', 'AT')
cd.AT = cd

cd = crim_data[grepl('^SK[0-9]{3}$', crim_data$geo), ]
cd = label_eurostat(cd, fix_duplicated = TRUE)
cd = subset(cd, unit == 'Number')
cd = subset(cd, freq == 'Annual')
cd = subset(cd, iccs == 'Intentional homicide')
cd = aggregate(list(values = cd$values), list(TIME_PERIOD = cd$TIME_PERIOD), sum)
mean.SK= mean(cd$values)
median.SK= median(cd$values)
colnames(cd) = c('year', 'SK')
cd.SK = cd

cd.AT$SK = cd.SK$SK
cd = cd.AT
cd

cat('AT mean, median: ', mean.AT, ',', median.AT, '\n')
cat('SK mean, median: ', mean.SK, ',', median.SK)

AT mean, median:  60.85714 , 59 
SK mean, median:  73.78571 , 75.5

Over 2008-2021, the median annual number of homicides is 59 in Austria and 75.5 in Slovakia. I used a one-sided Mann-Whitney U test, which does not assume normality, to test whether the Austrian counts are lower. Let \(\tilde{\mu}*\text{AT}\) and \(\tilde{\mu}*\text{SK}\) be the true respective median values.

\(H_0: \tilde{\mu}*\text{AT} = \tilde{\mu}*\text{SK}\)

\(H_A: \tilde{\mu}*\text{AT} < \tilde{\mu}*\text{SK}\)

wilcox.test(cd$AT, cd$SK, alternative='less', exact=F)

  Wilcoxon rank sum test with continuity correction

data:  cd$AT and cd$SK
W = 50, p-value = 0.01449
alternative hypothesis: true location shift is less than 0

The test rejects the null hypothesis at 5%: annual homicide counts are lower in Austria. Austria also has about 1.6 times the population of Slovakia, so the difference per capita is even larger.

hypothesis 3: normality of homicide data in Vienna

cd = crim_data
cd = subset(cd, unit == 'NR')
cd = subset(cd, freq == 'A')
cd = subset(cd, geo == 'AT13')
cd = label_eurostat(cd)
cd = subset(cd, iccs == 'Intentional homicide')

ggplot(cd, aes(x=values)) + xlab('Count of intentional homicide') + ylab('Density') +
    geom_histogram(aes(y=after_stat(density)), colour = 1, fill = 'white', bins=5) +
    stat_function(fun=dnorm, 
                  args=list(mean=mean(cd$values), sd=sd(cd$values)),
                  colour='red', lwd=2, linetype='dashed')

Finally, I checked whether the annual homicide counts in Vienna are consistent with a normal distribution, using a Shapiro-Wilk test and a Q-Q plot.

\(H_0\): The annual number of homicides in the Vienna region comes from a normal distribution.

\(H_A\): The annual number of homicides in the Vienna region does not come from a normal distribution.

shapiro.test(cd$values)

options(repr.plot.height=5)
ggplot(cd, aes(sample=values)) +
        stat_qq(distribution=qnorm, show.legend=T) +
        stat_qq_line(distribution=qnorm, show.legend=F)

  Shapiro-Wilk normality test

data:  cd$values
W = 0.94623, p-value = 0.5039

The test does not reject normality at 5%, and the points in the Q-Q plot stay close to the line. With this few observations the test has little power, so this is weak evidence for normality rather than a confirmation.