library(PxWebApiData) #Used to retrieve data from Statistics Norway
library(tidyverse)
library(janitor)
library(sf)
library(rvest)
library(plotly)
library(geojsonsf)
library(jsonlite)
library(readtext)
library(tinytable)STV2020 Final Project
1 Introduction
This project investigates social vulnerability across Oslo’s neighbourhoods by constructing a vulnerability index for each of the city’s 99 delbydeler (sub-districts). The aim is twofold: to identify which neighbourhoods can be characterised as the most vulnerable on a set of standardised socioeconomic indicators, and to examine how this pattern of vulnerability has changed between 2010 and 2023.
The project builds on Allvin, Gerell, and Skardhamar (2024)´s “Vulnerable Areas in Oslo – Toward ‘Swedish Conditions’?”, published in the Nordic Journal of Urban Studies. That study applies two established methodologies: the Danish threshold-based approach and the Swedish index-based approach. While the Danish method classifies areas as vulnerable when they meet at least two of four predefined criteria, the Swedish method produces a continuous vulnerability score by combining standardized indicators of social disadvantage.
In this paper, I adapt the Swedish method, originally developed by Guldåker, Hallin, Nilvall, and Gerell (2021). The method constructs a vulnerability index from three sub-indexes, each covering different dimensions of social vulnerability. The economic living-conditions sub-index combines education level, employment share, and income. The family-conditions sub-index combines the share of single-parent households and population density. The segregation sub-index is based on the share of the population with an immigrant background. Each indicator is standardised as a z-score relative to a reference population, capped between −1 and +3, and the three sub-indexes are averaged into a final score. Areas scoring more than two standard deviations above the mean are classified as vulnerable.
The central contribution of this paper is its temporal aspect. Whereas Allvin et al. (2024) analyze data from only the year of 2020, this study constructs the vulnerability index for each year between 2010 and 2023, making it possible to observe how the spatial distribution of vulnerability in Oslo has evolved over more than a decade.
2 Methodology
I follow the same three-dimensional structure as Guldåker et al. (2021), but with three methodological differences, each driven by data availability:
I collect data from 99 ‘delbydeler’ instead of 500+ ‘grunnkretser’, as data on such small geographic units is not available for the public.
I construct the reference population from all of Norway’s municipalities (kommuner), whereas Allvin et al. standardise against smaller units.
I use mean income rather than median income, as district-level median income was not available.
These differences will produce a result slightly skewed towards higher vulnerability indexes, compared to Allvin et al., as will be discussed in the conclusion.
The analysis draws on data from Oslo Kommunes Statistikkbank and Statistisk sentralbyrå (SSB). Most data is retrieved through API requests, and a helper function, fetch_data, is created at the beginning to handle these requests consistently. One dataset is obtained by scraping a saved HTML file, and one figure from Allvin et al.’s paper is extracted from a PDF using the readtext package.
Firstly, I build the national reference. I collect data from all of Norway’s municipalities over the time span 2010–2023 on six variables: population density, immigrant share, education level, mean gross income, employment share, and single-parent households. For each variable and each year, I compute a national mean and standard deviation, which Oslo districts later on is standardized against.
I then retrieve the same six variables for each of Oslo’s 99 delbydeler over the same time period. Joining the national reference to the district data produces a master dataset, from which a vulnerability score is calculated for every district in every year, following the Swedish method described above.
In the results section, I validate the indexes against the original study, comparing the average value of each variable within each vulnerability band across the two studies. I examine how vulnerability in Oslo as a whole has evolved over time, and I examine whether neighbourhoods have grown more or less equal in their levels of social disadvantage over time. Lastly, I create a heat map displaying each district and their vulnerability score, with an animation_slider, allowing us to see how it changes throughout the years.
3 Loading packages
4 Creating an API request function
As most data is retrieved through APIs, I start by creating a helper function that fetches and caches the data. The first time it runs, it queries the API, converts the result to a tibble, and saves it to disk; on every later call it loads that saved file instead, so the document doesn’t re-query the API each time it renders.
fetch_data <- function(url, filename, ...) {
if (file.exists(filename)) { # If we've already fetched and saved this data,
load(filename) # load it from disk instead of re-querying
return(data)
}
response <- ApiData(url, ...) # Otherwise query the API
data <- as_tibble(response[[1]]) # Keep the first frame as a tibble
save(data, file = filename) # Save so future renders can skip the API call
return(data)
}5 National data
5.1 Creating functions for the national data
As the national data is all retrieved from SSB and look relatively similar, I create a function for the shared cleaning steps the raw datasets needs. Due to Norway’s municipalities often changing names, values for a municipality during a year it didn’t exist will appear as zero. All rows with a value of zero is therefore filtered out. I then create a function that calculates the national mean and standard deviation per year.
# Data wrangling function
clean_national <- function(data, value_name) {
data %>%
filter(value != 0) %>% # Drop rows with a value of zero
rename(year = år, "{value_name}" := value) %>% # Rename the year column and the value column
mutate(year = as.numeric(year)) # Make year numeric
}
# Mean + SD function
national_stats <- function(data, column, name) {
data %>%
group_by(year) %>% # Compute the statistics separately for each year
summarise(
"{name}_national_mean" := mean({{ column }}, na.rm = TRUE),
"{name}_national_sd" := sd({{ column }}, na.rm = TRUE),
.groups = "drop") # Drop grouping so later operations aren't grouped by year
}In the following sections, national data on population density, immigrant population share, the share of the population with only a basic or no level of education, mean gross income, employment share and single-parent share will be retrieved using the fetch_data function, and the variable’s national mean and standard deviation for each year will be calculated using the clean_national and national_stats functions.
5.2 Population density
#Fetching the data from SSB
national_popdensity_raw <- fetch_data(
url = "https://data.ssb.no/api/pxwebapi/v2/tables/11342/data?lang=no&outputFormat=json-stat2&valuecodes[Region]=*&codelist[Region]=vs_Kommune&valuecodes[ContentsCode]=FolkeLandArealKm2&valuecodes[Tid]=2010,2011,2012,2013,2014,2015,2016,2017,2018,2019,2020,2021,2022,2023&heading=Tid,ContentsCode&stub=Region",
filename = "ssb_popdensity_national_meta.RData",
getDataByGET = TRUE,
Region = "Alle kommuner",
Statistikkvariabel = "Innbyggere per km^2 landareal",
år = as.character(2010:2023)
)#Calculate national mean + SD for each year
national_popdensity <- national_popdensity_raw %>%
clean_national(value_name = "pop_density") %>%
national_stats(column = pop_density, name = "popdensity")5.3 Immigrant population
#Fetching the data from SSB
national_immigrant_population_raw <- fetch_data(
url = "https://data.ssb.no/api/pxwebapi/v2/tables/09817/data?lang=no&outputFormat=json-stat2&valuecodes[Tid]=2010,2011,2012,2013,2014,2015,2016,2017,2018,2019,2020,2021,2022,2023&valuecodes[Region]=*&codelist[Region]=vs_Kommune&valuecodes[Landbakgrunn]=999&valuecodes[ContentsCode]=AndelBefolkning&valuecodes[InnvandrKat]=B-C&heading=Tid,ContentsCode,Landbakgrunn&stub=InnvandrKat,Region",
filename = "ssb_immigrant_national_meta.RData",
getDataByGET = TRUE,
Region = "Alle kommuner",
Statistikkvariabel = "Andel av befolkningen (prosent)",
Innvandringskategori = "Innvandrere og norskfødte med innvandrerforeldre",
Landbakgrunn = "Alle land",
år = as.character(2010:2023)
)#Calculate national mean + SD for each year
national_immigrant_population <- national_immigrant_population_raw %>%
clean_national(value_name = "immigrant_population_share") %>%
mutate(immigrant_population_share = immigrant_population_share / 100) %>%
national_stats(column = immigrant_population_share, name = "immigrant")5.4 Education level
#Fetching the data from SSB
national_education_raw <- fetch_data(
url = "https://data.ssb.no/api/pxwebapi/v2/tables/09429/data?lang=no&outputFormat=json-stat2&valuecodes[Tid]=2010,2011,2012,2013,2014,2015,2016,2017,2018,2019,2020,2021,2022,2023&valuecodes[Region]=*&codelist[Region]=vs_Kommune&valuecodes[Kjonn]=0&valuecodes[Nivaa]=01,09a&valuecodes[ContentsCode]=Personer&heading=ContentsCode,Tid,Kjonn&stub=Region,Nivaa",
filename = "ssb_education_national_meta.RData",
getDataByGET = TRUE,
bosted = TRUE,
Statistikkvariabel = "Personer 16 år og over",
kjønn = "Begge kjønn",
år = as.character(2010:2023)
)
#Data wrangling
national_education <- national_education_raw %>%
group_by(region, år) %>%
summarise(value = sum(value, na.rm = TRUE), .groups = "drop") %>%
clean_national(value_name = "basic_education")5.4.1 Population above 16 for reference
Data on the percentage of the population with only a basic education level (grunnskole) or no education don’t exist, it only comes in raw numbers. I therefore fetch another dataset on the population in each municipality, and filter out ages below 16.
national_population_raw <- fetch_data(
url = "https://data.ssb.no/api/pxwebapi/v2/tables/07459/data?lang=no&outputFormat=json-stat2&valuecodes[Tid]=2010,2011,2012,2013,2014,2015,2016,2017,2018,2019,2020,2021,2022,2023&valuecodes[Region]=*&codelist[Region]=vs_Kommun&valuecodes[Alder]=F303,F304&codelist[Alder]=agg_Funksjonell1a&valuecodes[Kjonn]=*&valuecodes[ContentsCode]=*&heading=ContentsCode,Tid&stub=Region,Alder,Kjonn",
filename = "ssb_pop_national_meta.RData",
getDataByGET = TRUE,
Region = "Alle kommuner",
Alder = "16 - 66 år", "67 år eller eldre",
Statistikkvariabel = "Personer",
år = as.character(2010:2023)
)
#Data wrangling
national_population <- national_population_raw %>%
group_by(region, år) %>%
summarise(value = sum(value, na.rm = TRUE), .groups = "drop") %>%
clean_national(value_name = "population_above_16")I add the population column to the education dataset, and calculate the share. I then calculate the national mean and standard deviation as usual.
#Calculate national mean + SD for each year
national_education <- national_education %>%
left_join(national_population) %>%
mutate(basic_education_share = basic_education / population_above_16) %>%
select(-population_above_16, -basic_education) %>%
national_stats(column = basic_education_share, name = "education")5.5 Mean gross income
#Fetching the data from SSB
national_income_raw <- fetch_data(
url = "https://data.ssb.no/api/pxwebapi/v2/tables/03068/data?lang=no&outputFormat=json-stat2&valuecodes[ContentsCode]=Bruttoinnt&valuecodes[Tid]=2010,2011,2012,2013,2014,2015,2016,2017,2018,2019,2020,2021,2022,2023&valuecodes[Region]=*&codelist[Region]=vs_Kommune&valuecodes[Kjonn]=0&heading=ContentsCode,Tid&stub=Region",
filename = "ssb_income_national_meta.RData",
getDataByGET = TRUE,
Region = "Alle kommuner",
Statistikkvariabel = "Gjennomsnittlig bruttoinntekt (kr)",
Kjønn = "Begge kjønn",
år = as.character(2010:2023)
)The mean gross income needs to be adjusted for inflation. I retrieve the monthly Consumer Price Index (CPI) from Statistics Norway for 2010–2023 and average the monthly values into an annual index for each year. I then compute an inflation factor for every year relative to the 2023 price level, dividing the 2023 index by each year’s index. Multiplying a nominal income by this factor expresses it in constant 2023 kroner, making incomes comparable across years.
#Adjusting for inflation
kpi <- fetch_data(
url = "https://data.ssb.no/api/pxwebapi/v2/tables/03013/data?lang=no&outputFormat=json-stat2&valuecodes[Tid]=*&valuecodes[ContentsCode]=KpiIndMnd&valuecodes[Konsumgrp]=*&codelist[Konsumgrp]=vs_CoiCop2016niva1&heading=Tid,Konsumgrp&stub=ContentsCode",
filename = "ssb_kpi.RData",
getDataByGET = TRUE,
år = as.character(2010:2023)
) %>%
rename(month = måned, kpi = value) %>%
mutate(year = as.numeric(str_sub(month, 1, 4)),
kpi = as.numeric(kpi)) %>%
group_by(year) %>%
summarise(kpi = mean(kpi, na.rm = TRUE), .groups = "drop") %>%
mutate(inflation_factor = kpi[year == 2023] / kpi)#Calculate national mean + SD for each year
national_income <- national_income_raw %>%
clean_national(value_name = "mean_income") %>%
left_join(select(kpi, year, inflation_factor)) %>%
mutate(mean_income_adjusted = mean_income * inflation_factor) %>%
national_stats(column = mean_income_adjusted, name = "income")5.6 Employment
#Fetching the data from SSB
national_employment_raw <- fetch_data(
url = "https://data.ssb.no/api/pxwebapi/v2/tables/06445/data?lang=no&outputFormat=json-stat2&valuecodes[Region]=*&codelist[Region]=vs_Kommune&valuecodes[Alder]=15-74&valuecodes[Kjonn]=0&valuecodes[Tid]=2010,2011,2012,2013,2014,2015,2016,2017,2018,2019,2020,2021,2022,2023&valuecodes[ContentsCode]=*&heading=ContentsCode,Tid,Alder&stub=Region,Kjonn",
filename = "ssb_employment_national_meta.RData",
getDataByGET = TRUE,
Region = "Alle kommuner",
Statistikkvariabel = "Sysselsatte i prosent av befolkningen 15-74 år",
Kjønn = "Begge kjønn",
år = as.character(2010:2023)
)#Calculate national mean + SD for each year
national_employment <- national_employment_raw %>%
clean_national(value_name = "employment_share") %>%
mutate(employment_share = employment_share / 100) %>%
national_stats(column = employment_share, name = "employment")5.8 National reference master dataset
Before moving on to the Oslo district data, I combine the six national baseline datasets into a single reference table using left_join(). Because each dataset shares only the year column as a common variable, left_join() automatically identifies year as the key to join on, so we do not need to specify it explicitly.
# Creating national reference table
national_reference <- national_income %>%
left_join(national_education) %>%
left_join(national_immigrant_population) %>%
left_join(national_employment) %>%
left_join(national_singleparent) %>%
left_join(national_popdensity)tt(national_reference)| year | income_national_mean | income_national_sd | education_national_mean | education_national_sd | immigrant_national_mean | immigrant_national_sd | employment_national_mean | employment_national_sd | singleparent_national_mean | singleparent_national_sd | popdensity_national_mean | popdensity_national_sd |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 2010 | 459020.8 | 45340.41 | 0.3542071 | 0.06275147 | 0.06668302 | 0.03344690 | 0.6920023 | 0.04225858 | 0.1848571 | 0.04258122 | 54.86885 | 149.3569 |
| 2011 | 478438.8 | 55235.43 | 0.3500420 | 0.06320357 | 0.07391977 | 0.03493844 | 0.6908628 | 0.04119981 | 0.1894438 | 0.04245465 | 55.64403 | 151.8303 |
| 2012 | 493378.7 | 47927.40 | 0.3472325 | 0.06263327 | 0.08217669 | 0.03663863 | 0.6863061 | 0.04032894 | 0.1906432 | 0.04071762 | 56.55869 | 154.4730 |
| 2013 | 501161.5 | 50440.79 | 0.3445251 | 0.06244816 | 0.09089252 | 0.03926621 | 0.6836776 | 0.04247103 | 0.1908964 | 0.04209645 | 57.48000 | 157.1259 |
| 2014 | 508534.4 | 51768.10 | 0.3204841 | 0.05807752 | 0.09761449 | 0.04090083 | 0.6812921 | 0.04203869 | 0.1841095 | 0.04131242 | 58.18824 | 159.2486 |
| 2015 | 519626.1 | 52807.80 | 0.3161619 | 0.05720365 | 0.10419509 | 0.04206023 | 0.6620327 | 0.04142315 | 0.1818178 | 0.04077416 | 58.88235 | 161.4410 |
| 2016 | 506168.9 | 46035.11 | 0.3118155 | 0.05612725 | 0.11117243 | 0.04342377 | 0.6598645 | 0.03951427 | 0.1802907 | 0.03901645 | 59.53412 | 163.2171 |
| 2017 | 510733.9 | 44112.35 | 0.3069777 | 0.05504893 | 0.11820188 | 0.04385693 | 0.6630258 | 0.03999908 | 0.1784423 | 0.03665064 | 59.52719 | 164.2942 |
| 2018 | 514406.3 | 43647.74 | 0.3004984 | 0.05433763 | 0.12301209 | 0.04458399 | 0.6669953 | 0.03926474 | 0.1780832 | 0.03686514 | 59.72315 | 165.9069 |
| 2019 | 525033.7 | 44136.24 | 0.2941135 | 0.05312750 | 0.12538531 | 0.04614858 | 0.6700379 | 0.03914676 | 0.1761308 | 0.03693385 | 60.14320 | 167.3235 |
| 2020 | 529575.0 | 44663.26 | 0.2886821 | 0.05211387 | 0.12637051 | 0.04788531 | 0.6638848 | 0.03793673 | 0.1797419 | 0.03659505 | 50.75637 | 131.2899 |
| 2021 | 556803.9 | 58614.47 | 0.2822689 | 0.05113183 | 0.12708764 | 0.05015384 | 0.6767978 | 0.03734459 | 0.1806005 | 0.03733805 | 51.02550 | 132.2353 |
| 2022 | 534983.2 | 44669.76 | 0.2772453 | 0.05002956 | 0.12906798 | 0.05229878 | 0.6784551 | 0.03741440 | 0.1783861 | 0.03803825 | 51.71795 | 133.6871 |
| 2023 | 537853.9 | 44707.24 | 0.2722432 | 0.04924478 | 0.14147921 | 0.05386170 | 0.6747500 | 0.03747678 | 0.1795109 | 0.03815258 | 52.46439 | 135.6937 |
6 Oslo district data
In the following section, I collect data from all of Oslo’s 99 sub-districts from the years 2010-2023 on the following variables: population density, immigrant population share, the share of the population with only a basic or no level of education, mean gross income, employment share and single-parent share. I use the fetch_data function here as well, but collect data from Oslo kommunes Statistikkbank.
6.1 Population of Oslo
#Fetching the data from Statistikkbanken
population_raw <- fetch_data(
url = "https://statistikkbanken.oslo.kommune.no/statbank/api/v1/no/db1/Befolkning/Folkemengde/OK-BEF002.px",
filename = "ssb_population_meta.RData",
bosted = TRUE,
aldersgruppe = TRUE ,
kjønn = TRUE ,
år = as.character(c(2010:2026))
)#Data wrangling
population <- population_raw %>%
group_by(bosted, år) %>%
summarise(
population = sum(value, na.rm = TRUE), # All ages
population_above_16 = sum(value[aldersgruppe %in% c("16-19 år", "20-29 år", "30-39 år", "40-49 år", "50-59 år", "60-66 år", "67-79 år", "80-89 år", "90 år +")], na.rm = TRUE), .groups = "drop") %>%
rename(district = bosted, year = år) %>%
mutate(year = as.numeric(year))6.2 Population density
To calculate the population density I dowload a GeoJSON file from Oslo Kommune’s website including all of Oslo’s sub-districts. I then calculate the area of each sub-district in km^2.
# Load GeoJSON file, transform to Norwegian coordinate system (metres) and calculate area in km^2
oslo_sf <- st_read("Data/Delbydeler_7090208274308320643.geojson", quiet = TRUE) %>%
st_transform(crs = 25833) %>% # Norwegian CRS, metres
mutate(area_km2 = as.numeric(st_area(geometry)) / 1000000) %>%
select(-FID, -kommunenum, -BYDEL, -BYDELSNAVN, -DELBYDEL) %>%
rename(district = DELBYDELSN) %>%
group_by(district) %>%
summarise(area_km2 = sum(area_km2),
geometry = st_union(geometry), .groups = "drop")I then join the area column and the population column together in a new dataset and calculate the population density for each sub-district.
# Calculating population density
population_density <- population %>%
select(district, year, population) %>%
left_join(
oslo_sf %>%
st_drop_geometry() %>%
select(district, area_km2)) %>% # Only bring in these two columns
mutate(pop_density = population / area_km2)6.3 Immigrant population
#Fetching data from Statistikkbanken
immigrant_population_raw <- fetch_data(
url = "https://statistikkbanken.oslo.kommune.no/statbank/api/v1/no/db1/Befolkning/Innvandrere/OK-BEF024.px",
filename = "ssb_befolkning_meta.RData",
bosted = TRUE,
alder = "9" ,
innvandringskategori = TRUE ,
år = as.character(c(2010:2023))
)immigrant_population <- immigrant_population_raw %>%
select(-EliminatedContents, -alder) %>%
rename(year = år, district = bosted, immigrant_population = value) %>%
mutate(year = as.numeric(year)) %>%
filter(innvandringskategori %in%
c("Innvandrer", "Norskfødt med innvandrerforeldre")) %>% # Keeping only the value for immigrants and children of immigrants
group_by(district, year) %>%
summarise(immigrant_population = sum(immigrant_population, na.rm = TRUE),
.groups = "drop") %>%
left_join(select(population, district, year, population)) %>% # Adding a population column
mutate(immigrant_share = immigrant_population / population) # Calculating the share of immigrantsThe immigrant share for the Sentrum district was not included in this dataset. Therefore, the data for Sentrum had to be retrieved by web scraping instead. I dowloaded and saved the HTML locally, so the analysis parses the saved file rather than the live page. This keeps the project reproducible: a scraped page may change or disappear over time, whereas a local snapshot always yields the same result.
The saved file is read with read_html(), which parses the HTML into a document object that can be navigated with selectors. The relevant table is then selected with html_element() using a CSS attribute selector, table[id^='12610'], which targets the element whose id begins with the Statistics Norway table number 12610.
Selecting by this identifier is more robust than relying on the table’s position on the page. Because html_element() (singular) returns a single node rather than a set, the matched table is passed directly to html_table(), which converts the HTML table into a data frame. The relevant rows and columns of this data frame are then extracted and reshaped into a tidy tibble with one row per year.
## Import HTML with immigrant data for Sentrum
# Parse the saved local HTML (kept local for reproducibility) and select the Sentrum table directly by its SSB table-number id prefix (12610).
page <- read_html("Data/12610_20260511-222724.html")
sentrum_raw <- page %>%
html_element(css = "table[id^='12610']") %>%
html_table()
# Build the Sentrum tibble
sentrum <- tibble(
district = "Sentrum",
year = as.integer(as.character(sentrum_raw[1, 2:15])),
immigrant_share = as.numeric(str_replace(as.character(sentrum_raw[4, 2:15]), ",", ".")) / 100
) %>%
arrange(year)# Adding the data to the immigrant_population dataset
immigrant_population <- immigrant_population %>%
filter(district != "Sentrum, Marka og uten registrert adresse") %>%
bind_rows(
sentrum %>%
mutate(
year = as.numeric(year),
immigrant_population = NA_real_,
population = NA_real_
)
)6.4 Education level
#Fetching the data from Statistikkbanken
education_raw <- fetch_data(
url = "https://statistikkbanken.oslo.kommune.no/statbank/api/v1/no/db1/Barnehage,%20skole%20og%20utdanning/Utdanningsniv%C3%A5/OK-UTD027.px",
filename = "ssb_utdanning_meta.RData",
bosted = TRUE,
aldersgruppe = TRUE,
kjønn = TRUE,
utdanningsnivå = c("Grunnskole", "Ingen utdanning/Uoppgitt utdanning"),
statistikkvariabel = TRUE,
år = as.character(c(2010:2023))
) %>%
mutate(år = as.numeric(as.character(år)))#Data wrangling
education <- education_raw %>%
select(-EliminatedContents, -statistikkvariabel, -NAstatus) %>%
rename(year = år, district = bosted, education_level = utdanningsnivå) %>% # Removes unnecessary columns and renames remaining columns
group_by(year, district) %>%
summarise(basic_education = sum(value, na.rm = TRUE), .groups = "drop")
#Calculating the share
education <- education %>%
left_join(select(population, district, year, population_above_16)) %>%
mutate(basic_education_share = basic_education / population_above_16) %>%
select(-basic_education, -population_above_16)6.5 Mean gross income
#Fetching data from Statistikkbanken
mean_income_raw <- fetch_data(
url = "https://statistikkbanken.oslo.kommune.no/statbank/api/v1/no/db1/Inntekt/OK-INN001.px",
filename = "ssb_inntekt_meta.RData",
Geografi = TRUE,
Alder = "Alder i alt",
Kjønn = "Begge kjønn",
Statistikkvariabel = "bruttoinntekt",
År = as.character(c(2010:2023))
)## Data wrangling
mean_income <- mean_income_raw %>%
select(-EliminatedContents, -Kjønn, -Alder, -Statistikkvariabel) %>%
rename(year = År, district = Geografi, mean_income = value) %>%
mutate(year = as.numeric(year)) %>%
filter(str_detect(district, "^\\d")) %>% # Using a regex pattern to filter the dataset so it keeps only rows where the district name starts with a digit, in this way removing "bydeler"
mutate(district = str_remove(district, "^\\d+ ")) %>% # Using a regex pattern to remove the numeric code and the following space from the start of each district name
left_join(select(kpi, year, inflation_factor)) %>%
mutate(mean_income = mean_income * inflation_factor) %>% # Adjusting for inflation
select(-inflation_factor)6.5.1 Building a district name lookup
Because the mean income data consists of both the district code and the district name, while a later dataset carries only the code, I build a lookup table that connects codes to names. I do this with stringr functions and regular expression patterns. str_extract() pulls out the leading digits as the code, str_remove() strips those digits to isolate the name, and str_detect() filters the result down to sub-districts by keeping only entries that begin with a digit.
district_lookup <- mean_income_raw %>%
mutate(
district_code = str_extract(Geografi, "^\\d+"), # Extracts one or more digits at the beginning of the string
district_name = str_remove(Geografi, "^\\d+ ") # Removes one or more digits at the start of the string, followed by a single space
) %>%
select(district_code, district_name, district_full = Geografi) %>%
distinct()
# Filter the dataset so it keeps only rows where the district name starts with a digit, in this way removing "bydeler" and keeping only "delbydeler"
district_lookup <- district_lookup %>%
filter(str_detect(district_full, "^\\d"))6.6 Employment
#Fetching data from Statistikkbanken
employment_raw <- fetch_data(
url = "https://statistikkbanken.oslo.kommune.no/statbank/api/v1/no/db1/Sysselsetting/OK-SYS002.px",
filename = "ssb_employment_meta.RData",
geografi = TRUE,
alder = "Alder i alt",
kjønn = "Begge kjønn",
statistikkvariabel = "andel sysselsatte",
år = as.character(c(2010:2023))
)## Data wrangling
employment <- employment_raw %>%
select(-EliminatedContents, -alder, -statistikkvariabel, -kjønn) %>%
rename(year = år,
district = geografi,
employment_share = value) %>%
mutate(employment_share = employment_share / 100,
year = as.numeric(year)) %>%
distinct()6.8 Distribution of variables
Before the data analysis, we can take a look at the distribution of the six variables used in the study. The box-plots show the minimum, maximum, median, the lower and upper quartile, as well as the outliers.
immigrant_population %>%
select(district, year, immigrant_share) %>%
left_join(select(education, district, year, basic_education_share)) %>%
left_join(select(mean_income, district, year, mean_income)) %>%
left_join(select(employment, district, year, employment_share)) %>%
left_join(select(single_parent_share, district, year, single_parent_share)) %>%
left_join(select(population_density, district, year, pop_density)) %>%
filter(year == 2023) %>% # single year: distribution only
pivot_longer(-c(district, year), names_to = "variable", values_to = "value") %>%
mutate(variable = recode(variable,
immigrant_share = "Immigrant share",
basic_education_share = "Basic education share",
employment_share = "Employment share",
single_parent_share = "Single-parent share",
mean_income = "Mean income (NOK)",
pop_density = "Population density"),
variable = factor(variable, levels = c( # set the facet order
"Population density",
"Immigrant share",
"Basic education share",
"Mean income (NOK)",
"Employment share",
"Single-parent share"))) %>%
ggplot(aes(y = value)) +
geom_boxplot(fill = "steelblue") +
facet_wrap(~ variable, scales = "free") + # free scales: variables differ in range
theme_bw() +
labs(title = "Distribution of each variable across Oslo districts, 2023") +
theme(axis.text.x = element_blank(), axis.ticks.x = element_blank(), axis.title.y = element_blank())
7 Data analysis
7.1 Master dataset
# Creating a master dataset with all six variables
master_data <- immigrant_population %>%
select(district, year, immigrant_share) %>%
left_join(education %>% select(district, year, basic_education_share)) %>%
left_join(mean_income %>% select(district, year, mean_income)) %>%
left_join(employment %>% select(district, year, employment_share)) %>%
left_join(single_parent_share %>% select(district, year, single_parent_share)) %>%
left_join(population_density %>% select(district, year, pop_density))7.2 Calculating vulnerability indexes
The vulnerability indexes are calculated by calculating a z-score: variables value - the national mean / the national standard deviation. The z-scores are capped between -1 and 3. The variables are then divided into three sub-indexes: economic living-conditions, family conditions and segregation conditions, which are then averaged. The z-scores for income and employment are flipped, so that higher scores mean more disadvantage.
master_data <- master_data %>%
left_join(national_reference)
master_data <- master_data %>%
mutate(
z_income = (mean_income - income_national_mean) / income_national_sd,
z_education = (basic_education_share - education_national_mean) / education_national_sd,
z_immigrant = (immigrant_share - immigrant_national_mean) / immigrant_national_sd,
z_employment = (employment_share - employment_national_mean) / employment_national_sd,
z_single_parent = (single_parent_share - singleparent_national_mean) / singleparent_national_sd,
z_density = (pop_density - popdensity_national_mean) / popdensity_national_sd
)
# Flip income and employment so that higher z-scores mean more disadvantage (low income and low employment indicate higher vulnerability)
master_data <- master_data %>%
mutate(
z_income = -1 * z_income,
z_employment = -1 * z_employment
)
# Cap all six z-scores to the [-1, 3] range used by the Swedish method.
master_data <- master_data %>%
mutate(across(starts_with("z_"), ~ pmax(-1, pmin(3, .x))))
# Build sub-indexes and final index
master_data <- master_data %>%
mutate(
subindex_economic = (z_income + z_education + z_employment) / 3,
subindex_family = (z_single_parent + z_density) / 2,
subindex_segregation = z_immigrant,
vulnerability_index = (subindex_economic + subindex_family + subindex_segregation) / 3
)
# Apply threshold
master_data <- master_data %>%
mutate(
vulnerability_category = case_when(
vulnerability_index < 0 ~ "Below average",
vulnerability_index < 1 ~ "Around average",
vulnerability_index < 2 ~ "At risk",
vulnerability_index >= 2 ~ "Vulnerable"
),
vulnerability_category = factor(
vulnerability_category,
levels = c("Below average", "Around average", "At risk", "Vulnerable")
)
)7.2.1 Distribution of vulnerability indexes
index_histogram <- master_data %>%
filter(year == 2023) %>%
ggplot(aes(x = vulnerability_index, fill = vulnerability_category)) +
geom_histogram(breaks = seq(floor(min(master_data$vulnerability_index, na.rm = TRUE)),
ceiling(max(master_data$vulnerability_index, na.rm = TRUE)),
by = 0.25),
closed = "left",
color = "grey40") +
scale_fill_manual(values = c(
"Below average" = "#FFFFFF",
"Around average" = "#B0B0B0",
"At risk" = "#606060",
"Vulnerable" = "#000000")) +
theme_bw() +
labs(x = "Vulnerability index",
y = "Number of districts",
fill = "Category")
ggplotly(index_histogram)7.3 Comparing to Allvin et al.’s study
In order to assess whether my adaptation of the Swedish method behaves as intended, I include a comparison between the 2020 band profile of my index and a corresponding table in Allvin et al.’s paper. As Allvin et al. only analyses data from 2020, constructing an equivalent table for that year allows a comparison, before extending the analysis across time. Allvin et al.’s table only include low education share, employment share and single-parent share, so I compare my own data to those three variables. The purpose of the comparison is to verify the gradient of disadvantage across the vulnerability bands; as areas move into higher bands, one would expect rising shares of low education, single-parent households and falling share of employment. However, some design differences prevent a numerical replication: the original study uses “grunnkretser” rather than “delbydeler” and standardizes against almost 14,000 national “grunnkretser”, rather than municipality-level statistics. Allvin et al. note, larger units of analysis tend to smooth out extreme values.
For my own data, I filter the master dataset to 2020, group the districts by vulnerability category, and compute the mean of each indicator within each band. To obtain the published figures, I read Allvin et al.’s paper from its PDF using readtext() and split the text into individual lines. I locate Table 2 by searching for its caption and the final data row, then slice out the lines in between. Because each value in the table is formatted as “mean (min–max)”, I use a regular expression to extract all decimal numbers from a row and keep every third one, which isolates the band means. Applying this to each indicator’s row produces a tidy table of the published values. Finally, I join my figures to the published ones by indicator and band, reshape the result so the four bands form columns and the two sources appear as rows, and present the comparison as a grouped table.
# This study's category means in 2020
mine <- master_data %>%
filter(year == 2020) %>% # The comparison year
group_by(category = vulnerability_category) %>% # Group by vulnerability band
summarise(
basic_education_share = mean(basic_education_share, na.rm = TRUE),
employment_share = mean(employment_share, na.rm = TRUE),
single_parent_share = mean(single_parent_share, na.rm = TRUE),
.groups = "drop"
) %>%
mutate(category = recode(as.character(category), # Relabel to the band names
"Below average" = "Index < 0",
"Around average" = "Index 0–1",
"At risk" = "Index 1–2",
"Vulnerable" = "Index > 2")) %>%
pivot_longer(-category, names_to = "indicator", values_to = "mine") # long form for the join#Reading the PDF file
study <- readtext(file = "Data/2024-vulnerable-areas-in-oslo-toward-swedish-conditions.pdf") # Reads each file and returns a data frame with one row per document
study_lines <- study$text %>% # Access the raw text string of the document
read_lines() # Splits on \n to produce a character vector where each element is one line
table_start <- study_lines %>%
str_detect(pattern = "Table 2. Average aggregated scores by Swedish classification") %>% # Locate Table 2
which()
table_end <- study_lines %>%
str_detect(pattern = "median_totinc") %>% # The income row (last data row)
which() # Its line number
table_lines <- study_lines[table_start:(table_end - 1)] # Retrieve only lines from the table
# Extraction function for one variable row
extract_means <- function(x) {
nums <- x %>%
str_extract_all(pattern = "\\d+\\.\\d+") %>% # All decimal numbers on the line
unlist() %>% # Flatten the list to one vector
as.numeric() # Convert the strings to numbers
nums[seq(1, length(nums), by = 3)] # Keep positions 1,4,7,10,13 (the means)
}
# The variables we compare, and the text identifying each table row
indicator_rows <- c(
basic_education_share = "Share with low",
employment_share = "Share employed",
single_parent_share = "Share single-parent"
)
# For each of the three indicators, we find the table line containing that indicator's label, extract its band means, and combine the results into one data frame of the values published by Allvin et al.
published <- lapply(indicator_rows, function(row_text) {
i <- table_lines %>%
str_detect(pattern = row_text) %>%
which()
means <- extract_means(table_lines[i]) # The five band means for that line
tibble(
category = c("Index < 0", "Index 0–1", "Index 1–2", "Index > 2"),
published = means[2:5] # Drop "Overall", keep the four bands
)
}) %>%
bind_rows(.id = "indicator") # List names become the `indicator` column# Join the two sources and reshape to the wide comparison layout
comparison <- published %>%
left_join(mine) %>% # Match on indicator + band
pivot_longer(c(published, mine), names_to = "source", values_to = "value") %>% # Stack the two sources
mutate(
source = recode(source, published = "Allvin et al.", mine = "This study"), # Readable source labels
indicator = factor(indicator,
levels = c("basic_education_share", "employment_share", "single_parent_share"),
labels = c("Low education share", "Employment share", "Single-parent share")),
category = factor(category, # Fix the band order
levels = c("Index < 0", "Index 0–1", "Index 1–2", "Index > 2")),
value = round(value, 3) # Round for display
) %>%
arrange(indicator, desc(source), category) %>% # "Allvin et al." above "This study"
pivot_wider(names_from = category, values_from = value) %>%
rename(" " = source) %>%
select(-indicator) # Bands become columns
tt(comparison,
width = 1) %>%
group_tt(i = list(
"Low education share" = 1,
"Employment share" = 3,
"Single-parent share" = 5
)) %>%
style_tt(i = c(1, 4, 7), bold = TRUE) | Index < 0 | Index 0–1 | Index 1–2 | Index > 2 | |
|---|---|---|---|---|
| Low education share | Low education share | Low education share | Low education share | Low education share |
| This study | 0.134 | 0.143 | 0.207 | 0.375 |
| Allvin et al. | 0.300 | 0.350 | 0.510 | 0.630 |
| Employment share | Employment share | Employment share | Employment share | Employment share |
| This study | 0.692 | 0.701 | 0.695 | 0.572 |
| Allvin et al. | 0.760 | 0.710 | 0.610 | 0.510 |
| Single-parent share | Single-parent share | Single-parent share | Single-parent share | Single-parent share |
| This study | 0.139 | 0.135 | 0.209 | 0.208 |
| Allvin et al. | 0.080 | 0.100 | 0.130 | 0.200 |
The comparison shows that two of the three indicators reproduce the gradient found in Allvin et al. Low education rises steadily across the bands (from 0.134 in the below-average band to 0.375 in the most vulnerable), mirroring the upward pattern in the original study, though at consistently lower absolute levels. Employment likewise declines toward the vulnerable band (from 0.692 to 0.572), consistent with lower employment signalling greater disadvantage. The single-parent indicator is less clear-cut: while it is highest in the most vulnerable band (0.208), the gradient is uneven, and the below-average band shows a higher share (0.139) than the slightly-elevated band (0.135), unlike the steadily rising pattern in Allvin et al.
Because this study classifies Oslo’s 99 delbydeler using an index standardised against Norwegian municipalities, while Allvin et al. standardise against several thousand smaller grunnkretser and the national distribution of those units, the band means therefore reflect different underlying sets of neighbourhoods, which is why they need not align even when each indicator is measured identically. Low education is the exception to this pattern. The education gap reflects a genuine difference in how the variable is measured, in which I have not been able to figure out.
7.4 Overall vulnerability over time
Before moving on to how vulnerability scores are distributed across Oslo, we can take a quick look at how vulnerability in Oslo as a whole has developed over the time period. Grouping the master dataset by year and taking the mean of the index collapses the 99 districts into a single value per year, and the line graph shows that Oslo as a whole has become less vulnerable between 2010 and 2023.
master_data %>%
group_by(year) %>%
summarise(mean_vulnerability = mean(vulnerability_index, na.rm = TRUE)) %>%
ggplot(aes(x = year, y = mean_vulnerability)) +
geom_line(linewidth = 1, color = "darkred") +
theme_bw() +
labs(title = "Average vulnerability in Oslo over time",
x = "Year", y = "Mean vulnerability index")
7.5 Inequality over time
However, the rising standard deviation indicates that Oslo’s districts have grown more unequal in their vulnerability over the period. The standard deviation measures how far districts spread around the average, so an increase means the gap between the least and most vulnerable districts has widened. Even if the city-wide average vulnerability decreased, the distribution around the mean has stretched, pointing to growing spatial inequality between neighbourhoods.
inequality_over_time <- master_data %>%
group_by(year) %>%
summarise(
sd_index = sd(vulnerability_index, na.rm = TRUE),
.groups = "drop")
inequality_over_time %>%
ggplot(aes(x = year, y = sd_index)) +
geom_line(linewidth = 1, colour = "darkred") +
theme_bw() +
labs(title = "Spread of vulnerability across Oslo districts, 2010-2023",
x = "Year", y = "Standard deviation of the index")
A boxplot summarises the distribution of a variable through five values: the median (the line inside the box), the first and third quartiles (the box edges, spanning the middle half of the observations), and the whiskers extending to the most extreme values within range. The height of the box therefore indicates how much the districts vary in a given year. I therefore create a boxplot for each year, displaying the distribution of the vulnerability indexes. The widening boxes across the years reinforce the finding that inequality between Oslo’s districts is increasing. Each box spans the middle half of the districts in a given year, so a taller box means the index values are more spread out; the districts are becoming less alike. This visual pattern corroborates the rising standard deviation reported above. Both describe the same underlying trend of growing spatial inequality, with the boxplot adding the detail of where in the distribution the spread is occurring.
# Distribution and development of the vulnerability index over time.
master_data %>%
ggplot(aes(x = factor(year), y = vulnerability_index)) +
geom_boxplot(fill = "grey80", outlier.size = 0.5) +
geom_hline(yintercept = 0, linetype = "dashed", colour = "grey40") +
geom_hline(yintercept = 2, linetype = "dashed", colour = "black") +
theme_bw() +
labs(title = "Distribution of the vulnerability index across Oslo, 2010–2023",
x = "Year", y = "Vulnerability index")
To see which districts drive the widening spread, I track the two ends of the distribution separately over time. For each year I calculate the 90th percentile of the vulnerability index (representing the more vulnerable districts) and the 10th percentile (the least vulnerable), and plot the two as separate lines. This is relevant because the standard deviation and boxplots show that inequality is increasing but not where the change comes from. A widening gap could reflect the most vulnerable districts getting worse, the least vulnerable improving, or both.
By following the top and bottom of the distribution side by side, this figure reveals that the widening gap is driven mainly by the least vulnerable districts rather than the most vulnerable. The 90th percentile (the vulnerable end) stays relatively steady across the period, while the 10th percentile (the least vulnerable end) declines. In other words, the increase in inequality reflects the least vulnerable districts moving toward lower vulnerability scores, rather than the most vulnerable districts becoming more vulnerable. The gap between Oslo’s neighbourhoods has widened less because the worst-off areas got worse and more because the best-off areas improved their relative position.
# The 90th percentile (the more vulnerable districts) and the 10th percentile (the least vulnerable districts) of the index over the years
master_data %>%
group_by(year) %>%
summarise(
p90 = quantile(vulnerability_index, 0.9, na.rm = TRUE),
p10 = quantile(vulnerability_index, 0.1, na.rm = TRUE),
.groups = "drop"
) %>%
pivot_longer(cols = c(p90, p10),
names_to = "percentile", values_to = "index") %>%
ggplot(aes(x = year, y = index, color = percentile)) +
geom_line(linewidth = 1) +
scale_color_manual(
values = c(p90 = "tomato2", p10 = "cornflowerblue"),
labels = c(p90 = "90th percentile (most vulnerable)",
p10 = "10th percentile (least vulnerable)")
) +
theme_bw() +
labs(title = "Most and least vulnerable Oslo districts, 2010-2023",
x = "Year", y = "Vulnerability index", color = NULL)
7.6 Heat map over time
Finally, I visualise the vulnerability index spatially as an animated choropleth map. I prepare the district polygons by dropping the forested Marka areas, reprojecting them to the coordinate system plotly expects, and converting them to GeoJSON, matching each polygon to its data through the district name. Each district is then shaded according to its vulnerability category, using a grayscale scheme that follows Allvin et al. (2024); white for districts below the national average, through two shades of grey, to black for the most vulnerable. In this way, the four bands appear as distinct blocks rather than a continuous gradient. A slider lets the map step through each year from 2010 to 2023, with one frame per year.
This spatial view is relevant because the index and the inequality measures show how much vulnerability and its spread have changed, but not where. The map adds the geographical dimension. It shows which parts of Oslo are most vulnerable, whether vulnerable districts cluster together, and how that spatial pattern shifts over time.
The map reveals a clear geographical pattern. Vulnerability is concentrated in the south and northeast parts of Oslo, with some slightly higher vulnerability also in the eastern parts of the city center. The western and northern districts remain consistently near the national average, or even below in the last few years. The most vulnerable districts cluster together in these northeastern and southern parts of the city rather than being scattered evenly, indicating that disadvantage in Oslo is spatially structured rather than randomly distributed. Allvin et al. conclude with the same findings of where vulnerability is concentrated, writing “These findings show an unsurprising overlap with the long-lasting patterns of residential segregation found in demographic research, showing concentrations of more resourceful residents in the western part of the city” (Alvin et al., 2024).
# Prepare the district polygons and convert them to GeoJSON. Drop the Marka areas, reproject to CRS 4326 (what plotly expects), turn the geometry into a GeoJSON string with sf_geojson(), and parse it into an editable R list with fromJSON(), the format plotly reads for the choropleth.
oslo_geojson <- oslo_sf %>%
filter(!str_detect(district, "Marka")) %>%
st_transform(crs = 4326) %>%
sf_geojson() %>%
fromJSON(simplifyVector = FALSE)
# For every district polygon, copy the district name from the feature's properties into its id field. Plotly uses this id to match each polygon on the map to its vulnerability value, linking the geometry to the data.
oslo_geojson$features <- lapply(oslo_geojson$features, function(f) {
f$id <- f$properties$district
f
})
# Build the data table with one row per district per year. st_drop_geometry() removes the polygons (they already live in oslo_geojson and are matched back through the id), left_join() adds the vulnerability scores, and mutate() reates a numeric category code and a hover label shown when the cursor is over a district.
oslo_map_data <- oslo_sf %>%
st_drop_geometry() %>% # Keep the data, drop the polygons
filter(!str_detect(district, "Marka")) %>% # Drop the Marka areas
left_join(
master_data %>% select(district, year, vulnerability_category, vulnerability_index)) %>% # Add the vulnerability scores
mutate(
cat_code = as.numeric(vulnerability_category),
hover = paste0(district, "\n", vulnerability_category,
"\nIndex: ", round(vulnerability_index, 2))
)
# Grayscale scheme following Allvin et al. (2024): white for districts below the average, through two shades of grey, to black for the most vulnerable.
category_levels <- c("Index < 0", "Index 0–1", "Index 1–2", "Index > 2")
category_colours <- c("#FFFFFF", "#B0B0B0", "#606060", "#000000")
# Build a stepped colorscale, so each of the four categories appears as one solid block rather than a gradient.
discrete_scale <- list(
list(0.00, "#FFFFFF"), list(0.25, "#FFFFFF"),
list(0.25, "#B0B0B0"), list(0.50, "#B0B0B0"),
list(0.50, "#606060"), list(0.75, "#606060"),
list(0.75, "#000000"), list(1.00, "#000000")
)
# Build the animated choropleth. add_trace() draws the map, geojson supplies the polygons, locations matches them to districts by id, z gives the value to colour by, and frame = year creates one animation frame per year.
oslo_graph <- plot_ly() %>%
add_trace(
type = "choroplethmapbox",
geojson = oslo_geojson, # The district polygons
locations = oslo_map_data$district, # Match polygons by district id
z = oslo_map_data$cat_code, # Value used for colour
frame = oslo_map_data$year, # One frame per year
zmin = 1, # Lowest category code
zmax = 4, # Highest category code
colorscale = discrete_scale, # The stepped white-to-black scale
text = oslo_map_data$hover, # Hover label text
hoverinfo = "text", # Show only the hover text
marker = list(line = list(width = 0.5, color = "grey70")), # District border lines
colorbar = list(
tickvals = c(1.375, 2.125, 2.875, 3.625), # Centre of each of the four bands
ticktext = category_levels, # Label each band
title = "Vulnerability",
len = 0.5,
thickness = 15
)
) %>%
layout(
title = "Vulnerability in Oslo, 2010–2023",
mapbox = list(style = "white-bg", zoom = 9.3, # Adds plain white background
center = list(lon = 10.75, lat = 59.9)) # Centre on Oslo
) %>%
config(displayModeBar = FALSE) %>% # Hide the plotly toolbar
animation_slider(currentvalue = list(prefix = "Year: ")) # Label the slider "Year: "
oslo_graph