library(tidyverse)
library(cbsodataR)

Doel en intro

Ga naar statline voor veiligheid en kerncijfers.

We willen het aantal verdachten per gemeente begrijpen vanuit de kerncijfers. Kijk in de kerncijfers welke cijfers eventueel relevant zijn voor het begrijpen van criminaliteit.

Optie: je mag ook een andere afhankelijke variabele kiezen, check in StatLine Nederland Regionaal welke statistieken per gemeente bekend zijn.

Gegevens: afhankelijke variabele

We halen eerst de gegevens op van verdachten per regio.

Noot: Om het sneller te maken filteren we meteen op periode, stedelijkheid, en kiezen we alleen gemeenten:

# Deze code krijgen studenten 'kado'
verdacht_meta = cbs_get_meta("81960NED")
verdacht = cbs_get_data("81960NED", Perioden="2020JJ00", StedelijkheidBuurt="T001086", RegioS=has_substring("GM"))

Selectie van relevant kolommen, verwijderen NAs, en koppelen metagegevens regios aan data

# In te vullen door studenten
regios = verdacht_meta$RegioS |> as_tibble() |> select(RegioS=Key, regio=Title)
verdacht = inner_join(verdacht, regios)
verdacht = select(verdacht, statcode=RegioS, gemeente=regio, verdachten=TotaalVerdachtenVanMisdrijven_8)
verdacht = na.omit(verdacht)

Gegevens: Koppelen met onafhankelijke variabele

Ophalen kerncijfers voor gemeentes in dezelfde periode:

# Kado
kerncijfers_meta = cbs_get_meta("70072ned")
kerncijfers = cbs_get_data("70072ned", RegioS=has_substring("GM"), Perioden="2020JJ00")

Selecteren van relevante kerncijfers (onafhankelijke variabele(n)) en koppelen aan verdachten:

dichtheid = kerncijfers |> select(
  statcode=RegioS, 
  populatie=TotaleBevolking_1, 
  dichtheid=Bevolkingsdichtheid_57, 
  inkomen=ParticuliereHuishoudensExclStudenten_122)
d = inner_join(verdacht, dichtheid)

Analyse: Exploratief

Maak een scatterplot (geom_point) van de afhankelijke variabele afgezet tegen de gekozen onafhankelijke variabele

ggplot(d, aes(x=dichtheid, y=verdachten)) + geom_point() +
  theme_bw() + ggtitle("Aantal verdachten per 1000 inwoners per gemeente")

Als bonus: probeer deze mooier te maken, met bv andere variabelen voor kleur en grootte en labels voor interessante gemeentes

ggplot(d, aes(x=dichtheid, y=verdachten, color=inkomen)) + geom_point(aes(size=populatie)) +
  scale_color_gradient2(trans="log", midpoint = log(median(d$inkomen, na.rm=T)), low = scales::muted("red"), mid=scales::muted("blue"), high="green") + 
  ggrepel::geom_text_repel(data=filter(d, verdachten>120 | dichtheid>4000 | inkomen > 80), aes(label=gemeente)) + 
  theme_bw() + ggtitle("Aantal verdachten per 1000 inwoners per gemeente")

Analyse: Visualisatie op een kaart

Om de gegevens op een kaart te plotten halen we eerst de kaart van Nederlandse gemeenten op:

library(sf)
library(gridExtra)
map = st_read("https://geodata.nationaalgeoregister.nl/cbsgebiedsindelingen/wfs?request=GetFeature&service=WFS&version=2.0.0&typeName=cbs_gemeente_2020_gegeneraliseerd&outputFormat=json")

## Reading layer `OGRGeoJSON' from data source 
##   `https://geodata.nationaalgeoregister.nl/cbsgebiedsindelingen/wfs?request=GetFeature&service=WFS&version=2.0.0&typeName=cbs_gemeente_2020_gegeneraliseerd&outputFormat=json' 
##   using driver `GeoJSON'
## Simple feature collection with 355 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 13565.4 ymin: 306846.2 xmax: 278026.1 ymax: 619352.4
## Projected CRS: Amersfoort / RD New

Deze koppelen we aan de gegevens per gemeente:

d2= inner_join(map, d)

Nu kunnen we een kaart plotten:

ggplot(d2) + geom_sf(aes(fill=verdachten))

Of ietsje mooier:

 ggplot(d2) + geom_sf(aes(fill=verdachten)) +
  theme_minimal() + 
  scale_fill_gradient(low="white", high=scales::muted("red"), name="Verdachten per 1000 inwonders") +
  ggtitle("Aantal verdachten")

ggplot(d2) + 
  geom_sf(aes(fill=dichtheid)) + 
  theme_minimal() + 
  scale_fill_gradient(low="white", high=scales::muted("blue"), trans="log", breaks=c(10, 100, 1000, 10000), name="Inwoners per km2") +
  ggtitle("Bevolkingsdichtheid")

Analyse: toetsend

Test de correlatie tussen de afhankelijke en onafhankelijke variabele

cor.test(d$verdachten, log(d$dichtheid))

## 
##  Pearson's product-moment correlation
## 
## data:  d$verdachten and log(d$dichtheid)
## t = 10.753, df = 352, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.4144934 0.5718415
## sample estimates:
##       cor 
## 0.4972452

Bonus: als je meerdere onafhankelijke variabelen had, test dan een lineair (regressie)model:

# Bonus: regressie
m = lm(data=d, verdachten ~ log(dichtheid) + log(inkomen) + log(populatie))
summary(m)

## 
## Call:
## lm(formula = verdachten ~ log(dichtheid) + log(inkomen) + log(populatie), 
##     data = d)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -36.315 -10.429  -1.076   8.137  51.265 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    207.2633    33.1204   6.258 1.16e-09 ***
## log(dichtheid)   7.8077     0.9437   8.274 2.87e-15 ***
## log(inkomen)   -73.7352     6.8799 -10.717  < 2e-16 ***
## log(populatie)   9.6191     1.3413   7.171 4.56e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.38 on 345 degrees of freedom
##   (5 observations deleted due to missingness)
## Multiple R-squared:  0.5718, Adjusted R-squared:  0.5681 
## F-statistic: 153.6 on 3 and 345 DF,  p-value: < 2.2e-16

# Bonus: plots van model
library(sjPlot)
tab_model(m)

	Totaal verdachten van misdrijven
Predictors	Estimates	CI	p
(Intercept)	207.26	142.12 – 272.41	<0.001
log(dichtheid)	7.81	5.95 – 9.66	<0.001
log(inkomen)	-73.74	-87.27 – -60.20	<0.001
log(populatie)	9.62	6.98 – 12.26	<0.001
Observations	349
R² / R² adjusted	0.572 / 0.568

plot_residuals(m)

Teamopdracht 1

2022-09-19