The significance of the region on the salary in Sweden, a comparison between different occupational groups

In my last post, I found that the region has a significant impact on the salary of engineers. Is the significance of the region unique to engineers or are there similar correlations in other occupational groups?

Statistics Sweden use NUTS (Nomenclature des Unités Territoriales Statistiques), which is the EU’s hierarchical regional division, to specify the regions.

The F-value from the Anova table is used as the single value to discriminate how much the region and salary correlates. For exploratory analysis, the Anova value seems good enough.

First, define libraries and functions.

library (tidyverse) 
## -- Attaching packages -------------------------------------- tidyverse 1.3.0 --
## v ggplot2 3.2.1     v purrr   0.3.3
## v tibble  2.1.3     v dplyr   0.8.3
## v tidyr   1.0.2     v stringr 1.4.0
## v readr   1.3.1     v forcats 0.4.0
## -- Conflicts ----------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library (broom) 
library (car)
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
## The following object is masked from 'package:purrr':
## 
##     some
library (swemaps) # devtools::install_github('reinholdsson/swemaps')
library(sjPlot)
## Registered S3 methods overwritten by 'lme4':
##   method                          from
##   cooks.distance.influence.merMod car 
##   influence.merMod                car 
##   dfbeta.influence.merMod         car 
##   dfbetas.influence.merMod        car
## Install package "strengejacke" from GitHub (`devtools::install_github("strengejacke/strengejacke")`) to load all sj-packages at once!
readfile <- function (file1){  
  read_csv (file1, col_types = cols(), locale = readr::locale (encoding = "latin1"), na = c("..", "NA")) %>%  
    gather (starts_with("19"), starts_with("20"), key = "year", value = salary) %>%  
    drop_na() %>%  
    mutate (year_n = parse_number (year))
}
nuts <- read.csv("nuts.csv") %>%
  mutate(NUTS2_sh = substr(NUTS2, 1, 4))
nuts %>% 
  distinct (NUTS2_en) %>%
  knitr::kable(
    booktabs = TRUE,
    caption = 'Nomenclature des Unités Territoriales Statistiques (NUTS)')

map_ln_n <- map_ln %>%
  mutate(lnkod_n = as.numeric(lnkod)) 

The data table is downloaded from Statistics Sweden. It is saved as a comma-delimited file without heading, 000000CG.csv, http://www.statistikdatabasen.scb.se/pxweb/en/ssd/.

The table: Average basic salary, monthly salary and women´s salary as a percentage of men´s salary by region, sector, occupational group (SSYK 2012) and sex. Year 2014 – 2018 Monthly salary All sectors

In the plot and tables, you can also find information on how the increase in salaries per year for each occupational group is affected when the interactions are taken into account.

tb <- readfile ("000000CG.csv") %>%
  left_join(nuts %>% distinct (NUTS2_en, NUTS2_sh), by = c("region" = "NUTS2_en")) 
## Warning: Column `region`/`NUTS2_en` joining character vector and factor,
## coercing into character vector
tb_map <- readfile ("000000CG.csv") %>%
  left_join(nuts, by = c("region" = "NUTS2_en")) %>%  
  right_join(map_ln_n, by = c("Länskod" = "lnkod_n"))
## Warning: Column `region`/`NUTS2_en` joining character vector and factor,
## coercing into character vector
summary_table = vector()
anova_table = vector()
for (i in unique(tb$`occuptional  (SSYK 2012)`)){
  temp <- filter(tb, `occuptional  (SSYK 2012)` == i)
  if (dim(temp)[1] > 75){
    model <- lm(log(salary) ~ region + sex + year_n, data = temp)
    summary_table <- rbind (summary_table, mutate (tidy (summary (model)), ssyk = i, interaction = "none"))
    anova_table <- rbind (anova_table, mutate (tidy (Anova (model, type = 2)), ssyk = i, interaction = "none"))  
  
    model <- lm(log(salary) ~ region * sex + year_n, data = temp)
    summary_table <- rbind (summary_table, mutate (tidy (summary (model)), ssyk = i, interaction = "region and sex"))
    anova_table <- rbind (anova_table, mutate (tidy (Anova (model, type = 2)), ssyk = i, interaction = "region and sex"))
    
    model <- lm(log(salary) ~ region * year_n + sex, data = temp)
    summary_table <- rbind (summary_table, mutate (tidy (summary (model)), ssyk = i, interaction = "region and year"))
    anova_table <- rbind (anova_table, mutate (tidy (Anova (model, type = 2)), ssyk = i, interaction = "region and year"))   
    
    model <- lm(log(salary) ~ region * year_n * sex, data = temp)
    summary_table <- rbind (summary_table, mutate (tidy (summary (model)), ssyk = i, interaction = "region, year and sex"))
    anova_table <- rbind (anova_table, mutate (tidy (Anova (model, type = 2)), ssyk = i, interaction = "region, year and sex"))     
  }
}
## Note: model has aliased coefficients
##       sums of squares computed by model comparison
## Note: model has aliased coefficients
##       sums of squares computed by model comparison
anova_table <- anova_table %>% rowwise() %>% mutate(contcol = str_count(term, ":")) 

summary_table <- summary_table %>% rowwise() %>% mutate(contcol = str_count(term, ":"))

merge(summary_table, anova_table, by = c("ssyk", "interaction"), all = TRUE) %>%
  filter (term.x == "year_n") %>%
  filter (term.y == "region") %>%    
  filter (interaction == "none") %>%
  
  mutate (estimate = (exp(estimate) - 1) * 100) %>%  
  ggplot () +
    geom_point (mapping = aes(x = estimate, y = statistic.y, colour = interaction)) +
    labs(
      x = "Increase in salaries (% / year)",
      y = "F-value for interaction"
    )   

Figure 1: The significance of the region on the salary in Sweden, a comparison between different occupational groups, Year 2014 – 2018

merge(summary_table, anova_table, by = c("ssyk", "interaction"), all = TRUE) %>%
  filter (term.x == "year_n") %>%
  filter (contcol.y > 0) %>%    
  # only look at the interactions between all three variables in the case with interaction region, year and sex
  filter (!(contcol.y == 1 & interaction == "region, year and sex")) %>% 
  
  mutate (estimate = (exp(estimate) - 1) * 100) %>%  
  ggplot () +
    geom_point (mapping = aes(x = estimate, y = statistic.y, colour = interaction)) +
    labs(
      x = "Increase in salaries (% / year)",
      y = "F-value for interaction"
    )   

Figure 2: The significance of the interaction between region, year and sex on the salary in Sweden, a comparison between different occupational groups, Year 2014 – 2018

The tables with all occupational groups sorted by F-value in descending order.

merge(summary_table, anova_table, c("ssyk", "interaction"), all = TRUE) %>%
  filter (term.x == "year_n") %>%
  filter (term.y == "region") %>%   
  filter (interaction == "none") %>%
  mutate (estimate = (exp(estimate) - 1) * 100) %>%  
  select (ssyk, estimate, statistic.y, interaction) %>%
  rename (`F-value for age` = statistic.y) %>%
  rename (`Increase in salary` = estimate) %>%
  arrange (desc (`F-value for age`)) %>%
  knitr::kable(
    booktabs = TRUE,
    caption = 'Correlation for F-value (region) and the yearly increase in salaries')

merge(summary_table, anova_table, c("ssyk", "interaction"), all = TRUE) %>%
  filter (term.x == "year_n") %>%
  filter (contcol.y > 0) %>%   
  filter (interaction == "region and sex") %>%
  mutate (estimate = (exp(estimate) - 1) * 100) %>%  
  select (ssyk, estimate, statistic.y, interaction) %>%
  rename (`F-value for age` = statistic.y) %>%
  rename (`Increase in salary` = estimate) %>%
  arrange (desc (`F-value for age`)) %>%
  knitr::kable(
    booktabs = TRUE,
    caption = 'Correlation for F-value (region and sex) and the yearly increase in salaries')

merge(summary_table, anova_table, c("ssyk", "interaction"), all = TRUE) %>%
  filter (term.x == "year_n") %>%
  filter (contcol.y > 0) %>%   
  filter (interaction == "region and year") %>%
  mutate (estimate = (exp(estimate) - 1) * 100) %>%  
  select (ssyk, estimate, statistic.y, interaction) %>%
  rename (`F-value for age` = statistic.y) %>%
  rename (`Increase in salary` = estimate) %>%
  arrange (desc (`F-value for age`)) %>%
  knitr::kable(
    booktabs = TRUE,
    caption = 'Correlation for F-value (region and year) and the yearly increase in salaries')

merge(summary_table, anova_table, c("ssyk", "interaction"), all = TRUE) %>%
  filter (term.x == "year_n") %>%
  filter (contcol.y > 1) %>%   
  filter (interaction == "region, year and sex") %>%
  filter (!(contcol.y == 1 & interaction == "region, year and sex")) %>%   
  mutate (estimate = (exp(estimate) - 1) * 100) %>%  
  select (ssyk, estimate, statistic.y, interaction) %>%
  rename (`F-value for age` = statistic.y) %>%
  rename (`Increase in salary` = estimate) %>%
  arrange (desc (`F-value for age`)) %>%
  knitr::kable(
    booktabs = TRUE,
    caption = 'Correlation for F-value (region, year and sex) and the yearly increase in salaries')

Let’s check what we have found.

temp <- tb %>%
  filter(`occuptional  (SSYK 2012)` == "234 Primary- and pre-school teachers")
 
model <-lm (log(salary) ~ year_n + sex + NUTS2_sh, data = temp) 
 
plot_model(model, type = "pred", terms = c("NUTS2_sh"))
## Model has log-transformed response. Back-transforming predictions to original response scale. Standard errors are still on the log-scale.

Figure 3: Highest F-value region, Primary- and pre-school teachers

tb_map %>%
    filter(`occuptional  (SSYK 2012)` == "234 Primary- and pre-school teachers") %>%
  ggplot() +
    geom_polygon(mapping = aes(x = ggplot_long, y = ggplot_lat, group = lnkod, fill = salary)) +
      facet_grid(. ~ year) + 
    coord_equal() 

Figure 4: Highest F-value region, Primary- and pre-school teachers

temp <- tb %>%
  filter(`occuptional  (SSYK 2012)` == "911 Cleaners and helpers")
 
model <-lm (log(salary) ~ year_n + sex + NUTS2_sh, data = temp) 
 
plot_model(model, type = "pred", terms = c("NUTS2_sh"))
## Model has log-transformed response. Back-transforming predictions to original response scale. Standard errors are still on the log-scale.

Figure 5: Lowest F-value region, Cleaners and helpers

temp <- tb %>%
  filter(`occuptional  (SSYK 2012)` == "331 Financial and accounting associate professionals")
 
model <-lm (log(salary) ~ year_n + sex * NUTS2_sh, data = temp) 
 
plot_model(model, type = "pred", terms = c("NUTS2_sh", "sex"))
## Model has log-transformed response. Back-transforming predictions to original response scale. Standard errors are still on the log-scale.

Figure 6: Highest F-value interaction gender and region, Financial and accounting associate professionals

tb_map %>%
    filter(`occuptional  (SSYK 2012)` == "331 Financial and accounting associate professionals") %>%
  filter (sex == "men") %>%
  ggplot() +
    geom_polygon(mapping = aes(x = ggplot_long, y = ggplot_lat, group = lnkod, fill = salary)) +
      facet_grid(. ~ year) + 
    coord_equal() 

Figure 7: Highest F-value interaction gender and region, Financial and accounting associate professionals

tb_map %>%
    filter(`occuptional  (SSYK 2012)` == "331 Financial and accounting associate professionals") %>%
  filter (sex == "women") %>%  
  ggplot() +
    geom_polygon(mapping = aes(x = ggplot_long, y = ggplot_lat, group = lnkod, fill = salary)) +
      facet_grid(. ~ year) + 
    coord_equal() 

Figure 8: Highest F-value interaction gender and region, Financial and accounting associate professionals

temp <- tb %>%
  filter(`occuptional  (SSYK 2012)` == "222 Nursing professionals")
 
model <-lm (log(salary) ~ year_n + sex * NUTS2_sh, data = temp) 
 
plot_model(model, type = "pred", terms = c("NUTS2_sh", "sex"))
## Model has log-transformed response. Back-transforming predictions to original response scale. Standard errors are still on the log-scale.

Figure 9: Lowest F-value interaction gender and region, Nursing professionals

temp <- tb %>%
  filter(`occuptional  (SSYK 2012)` == "222 Nursing professionals")
 
model <-lm (log(salary) ~ year_n * NUTS2_sh + sex , data = temp) 
 
plot_model(model, type = "pred", terms = c("NUTS2_sh", "year_n"))
## Model has log-transformed response. Back-transforming predictions to original response scale. Standard errors are still on the log-scale.

Figure 10: Highest F-value interaction year and region, Nursing professionals

tb_map %>%
  filter(`occuptional  (SSYK 2012)` == "222 Nursing professionals") %>%
  ggplot() +
    geom_polygon(mapping = aes(x = ggplot_long, y = ggplot_lat, group = lnkod, fill = salary)) +
      facet_grid(. ~ year) + 
    coord_equal() 

Figure 11: Highest F-value interaction year and region, Nursing professionals

temp <- tb %>%
  filter(`occuptional  (SSYK 2012)` == "541 Other surveillance and security workers")
 
model <-lm (log(salary) ~ year_n * NUTS2_sh + sex , data = temp) 
 
plot_model(model, type = "pred", terms = c("NUTS2_sh", "year_n"))
## Model has log-transformed response. Back-transforming predictions to original response scale. Standard errors are still on the log-scale.

Figure 12: Lowest F-value interaction year and region, Other surveillance and security workers

temp <- tb %>%
  filter(`occuptional  (SSYK 2012)` == "531 Child care workers and teachers aides")
 
model <-lm (log(salary) ~ year_n * NUTS2_sh * sex , data = temp) 
 
plot_model(model, type = "pred", terms = c("NUTS2_sh", "year_n", "sex"))
## Model has log-transformed response. Back-transforming predictions to original response scale. Standard errors are still on the log-scale.

Figure 13: Highest F-value interaction region, year and gender, Child care workers and teachers aides

tb_map %>%
  filter(`occuptional  (SSYK 2012)` == "531 Child care workers and teachers aides") %>%
  filter (sex == "men") %>%
  ggplot() +
    geom_polygon(mapping = aes(x = ggplot_long, y = ggplot_lat, group = lnkod, fill = salary)) +
      facet_grid(. ~ year) + 
    coord_equal() 

Figure 14: Highest F-value interaction region, year and gender, Child care workers and teachers aides

tb_map %>%
    filter(`occuptional  (SSYK 2012)` == "531 Child care workers and teachers aides") %>%
  filter (sex == "women") %>%  
  ggplot() +
    geom_polygon(mapping = aes(x = ggplot_long, y = ggplot_lat, group = lnkod, fill = salary)) +
      facet_grid(. ~ year) + 
    coord_equal() 

Figure 15: Highest F-value interaction region, year and gender, Child care workers and teachers aides

temp <- tb %>%
  filter(`occuptional  (SSYK 2012)` == "336 Police officers")
 
model <-lm (log(salary) ~ year_n * NUTS2_sh * sex , data = temp) 
 
plot_model(model, type = "pred", terms = c("NUTS2_sh", "year_n", "sex"))
## Model has log-transformed response. Back-transforming predictions to original response scale. Standard errors are still on the log-scale.

Figure 16: Lowest F-value interaction region, year and gender, Police officers