So far I have analysed the effect of experience, education, gender and year on the salary of engineers in Sweden. In this post, I will have a look at the effect of the region on the salary of engineers in Sweden.

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

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))

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

We expect that the region is an important factor in salaries. As a null hypothesis, we assume that the region is not related to the salary and examine if we can reject this hypothesis with the data from Statistics Sweden.

tb <- readfile ("000000CG.csv") %>%
  filter (`occuptional  (SSYK 2012)` == "214 Engineering professionals") %>% 
  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") %>%
  filter (`occuptional  (SSYK 2012)` == "214 Engineering professionals") %>%
  left_join(nuts, by = c("region" = "NUTS2_en")) 
## Warning: Column `region`/`NUTS2_en` joining character vector and factor,
## coercing into character vector
tb_map %>%
  filter (sex == "men") %>%
  right_join(map_ln_n, by = c("Länskod" = "lnkod_n")) %>%
  ggplot() +
    geom_polygon(mapping = aes(x = ggplot_long, y = ggplot_lat, group = lnkod, fill = salary)) +
      facet_grid(. ~ year) + 
    coord_equal() 

Figure 1: SSYK 214, Architects, engineers and related professionals, men, Year 2014 – 2018

tb_map %>%
  filter (sex == "women") %>%
  right_join(map_ln_n, by = c("Länskod" = "lnkod_n")) %>%
  ggplot() +
    geom_polygon(mapping = aes(x = ggplot_long, y = ggplot_lat, group = lnkod, fill = salary)) +
      facet_grid(. ~ year) + 
    coord_equal() 

Figure 2: SSYK 214, Architects, engineers and related professionals, women, Year 2014 – 2018

tb %>%
  ggplot () +  
    geom_point (mapping = aes(x = year_n, y = salary, colour = region, shape=sex)) + 
  labs(
    x = "Year",
    y = "Salary (SEK/month)"
  )

Figure 3: SSYK 214, Architects, engineers and related professionals, Year 2014 – 2018

The F-value from the Anova table for the region is 40 (Pr(>F) < 2.2e-16), sufficient for rejecting the null hypothesis that the region has no effect on the salary holding year as constant. The adjusted R-squared value is 0,882 implying a good fit of the model.

model <-lm (log(salary) ~ year_n + sex + region , data = tb)

tb <- bind_cols(tb, as_tibble(exp(predict(model, tb, interval = "confidence"))))

tb %>%
  ggplot () +  
    geom_point (mapping = aes(x = year_n,y = fit, colour = region, shape=sex)) + 
  labs(
    x = "Year",
    y = "Salary (SEK/month)"
  )

Figure 4: Model fit, SSYK 214, Architects, engineers and related professionals, Year 2014 – 2018

summary(model) %>%  
  tidy() %>%
  knitr::kable( 
  booktabs = TRUE,
  caption = 'Summary from linear model fit')

summary(model)$r.squared
## [1] 0.8817871
Anova(model, type=2) %>% 
  tidy() %>% 
  knitr::kable( 
  booktabs = TRUE,
  caption = 'Anova report from linear model fit')

Let’s check what we have found.

model <-lm (log(salary) ~ year_n + sex + NUTS2_sh , data = tb) 
 
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: SSYK 214, Architects, engineers and related professionals, Year 2014 - 2018

plot(model, which = 1)

Figure 6: SSYK 214, Architects, engineers and related professionals, Year 2014 - 2018

tb[38,]
## # A tibble: 1 x 11
##   region sector `occuptional  (~ sex   year  salary year_n NUTS2_sh    fit
##   <chr>  <chr>  <chr>            <chr> <chr>  <dbl>  <dbl> <chr>     <dbl>
## 1 SE21 ~ 0 all~ 214 Engineering~ women 2016   34700   2016 SE21     38698.
## # ... with 2 more variables: lwr <dbl>, upr <dbl>
tb[27,]
## # A tibble: 1 x 11
##   region sector `occuptional  (~ sex   year  salary year_n NUTS2_sh    fit
##   <chr>  <chr>  <chr>            <chr> <chr>  <dbl>  <dbl> <chr>     <dbl>
## 1 SE31 ~ 0 all~ 214 Engineering~ men   2015   38100   2015 SE31     41616.
## # ... with 2 more variables: lwr <dbl>, upr <dbl>
tb[5,]
## # A tibble: 1 x 11
##   region sector `occuptional  (~ sex   year  salary year_n NUTS2_sh    fit
##   <chr>  <chr>  <chr>            <chr> <chr>  <dbl>  <dbl> <chr>     <dbl>
## 1 SE21 ~ 0 all~ 214 Engineering~ men   2014   41500   2014 SE21     39442.
## # ... with 2 more variables: lwr <dbl>, upr <dbl>

Also examine the interaction between gender and region. The F-value from the Anova table for the interaction is 2,7 (Pr(>F) < 0.017) implying a relation to 95 % significance.

model <-lm (log(salary) ~ year_n + sex * NUTS2_sh , data = tb) 

summary(model) %>%  
  tidy() %>%
  knitr::kable( 
  booktabs = TRUE,
  caption = 'Summary from linear model fit')

summary(model)$r.squared
## [1] 0.908906
Anova(model, type=2) %>% 
  tidy() %>% 
  knitr::kable( 
  booktabs = TRUE,
  caption = 'Anova report from linear model fit')

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 7: SSYK 214, Architects, engineers and related professionals, Year 2014 - 2018

plot(model, which = 1)

Figure 8: SSYK 214, Architects, engineers and related professionals, Year 2014 - 2018

tb[37,]
## # A tibble: 1 x 11
##   region sector `occuptional  (~ sex   year  salary year_n NUTS2_sh    fit
##   <chr>  <chr>  <chr>            <chr> <chr>  <dbl>  <dbl> <chr>     <dbl>
## 1 SE21 ~ 0 all~ 214 Engineering~ men   2016   40100   2016 SE21     41037.
## # ... with 2 more variables: lwr <dbl>, upr <dbl>

And the interaction between year and region. The F-value from the Anova table for the region is 0,57 (Pr(>F) < 0,77) implying no significant relation.

model <-lm (log(salary) ~ year_n * NUTS2_sh + sex , data = tb) 
 
summary(model) %>%  
  tidy() %>%
  knitr::kable( 
  booktabs = TRUE,
  caption = 'Summary from linear model fit')

summary(model)$r.squared
## [1] 0.8888956
Anova(model, type=2) %>% 
  tidy() %>% 
  knitr::kable( 
  booktabs = TRUE,
  caption = 'Anova report from linear model fit')

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 9: SSYK 214, Architects, engineers and related professionals, Year 2014 - 2018

plot(model, which = 1)

Figure 10: SSYK 214, Architects, engineers and related professionals, Year 2014 - 2018

Finally the interaction between gender, year and region, the only significant interaction is between the region and gender, F-value from the Anova table for the interaction is 2,4 (Pr(>F) < 0.033)

model <-lm (log(salary) ~ year_n * NUTS2_sh * sex , data = tb) 
 
summary(model) %>%  
  tidy() %>%
  knitr::kable( 
  booktabs = TRUE,
  caption = 'Summary from linear model fit')

summary(model)$r.squared
## [1] 0.9231133
Anova(model, type=2) %>% 
  tidy() %>% 
  knitr::kable( 
  booktabs = TRUE,
  caption = 'Anova report from linear model fit')

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 11: SSYK 214, Architects, engineers and related professionals, Year 2014 - 2018

plot(model, which = 1)

Figure 12: SSYK 214, Architects, engineers and related professionals, Year 2014 - 2018