The significance of the sector on the salary in Sweden, a comparison between different occupational groups, part 2

In my last post, I examined the significance of the sector on the salary for different occupational groups using statistics from different regions. In previous posts I have shown a correlation between the salary and experience and also salary and education, In this post, I will examine the correlation between salary and sector using statistics for age.

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

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

I have renamed the file to 000000D2_sector.csv because the filename 000000D2.csv was used in a previous post.

The table: Average basic salary, monthly salary and women´s salary as a percentage of men´s salary by sector, occupational group (SSYK 2012), sex and age. Year 2014 – 2018 Monthly salary 1-3 public sector 4-5 private sector

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("000000D2_sector.csv") %>%
  rowwise() %>%
  mutate(age_l = unlist(lapply(strsplit(substr(age, 1, 5), "-"), strtoi))[1]) %>%
  rowwise() %>%
  mutate(age_h = unlist(lapply(strsplit(substr(age, 1, 5), "-"), strtoi))[2]) %>%
  mutate(age_n = (age_l + age_h) / 2)

summary_table = 0
anova_table = 0

for (i in unique(tb$`occuptional  (SSYK 2012)`)){
  temp <- filter(tb, `occuptional  (SSYK 2012)` == i)
  if (dim(temp)[1] > 90){
    model <- lm(log(salary) ~ poly(age_n, 3) + sex + year_n + sector, 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) ~ poly(age_n, 3) * sector + sex + year_n, data = temp)
    summary_table <- rbind (summary_table, mutate (tidy (summary (model)), ssyk = i, interaction = "sector and age"))
    anova_table <- rbind (anova_table, mutate (tidy (Anova (model, type = 2)), ssyk = i, interaction = "sector and age"))  
    
    model <- lm(log(salary) ~ poly(age_n, 3) + sector * sex + year_n, data = temp)
    summary_table <- rbind (summary_table, mutate (tidy (summary (model)), ssyk = i, interaction = "sector and sex"))
    anova_table <- rbind (anova_table, mutate (tidy (Anova (model, type = 2)), ssyk = i, interaction = "sector and sex"))  
    
    model <- lm(log(salary) ~ poly(age_n, 3) +  year_n * sector + sex, data = temp)
    summary_table <- rbind (summary_table, mutate (tidy (summary (model)), ssyk = i, interaction = "sector and year"))
    anova_table <- rbind (anova_table, mutate (tidy (Anova (model, type = 2)), ssyk = i, interaction = "sector and year"))  
    
    model <- lm(log(salary) ~ poly(age_n, 3) * sector * sex * year_n, data = temp)
    summary_table <- rbind (summary_table, mutate (tidy (summary (model)), ssyk = i, interaction = "sector, year, age and sex"))
    anova_table <- rbind (anova_table, mutate (tidy (Anova (model, type = 2)), ssyk = i, interaction = "sector, year, age 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
## Note: model has aliased coefficients
##       sums of squares computed by model comparison
## 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 == "sector") %>%    
  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 sector"
    )   

Figure 1: The significance of the sector 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 four variables in the case with interaction sector, year, age and sex
  filter (!(contcol.y < 3 & interaction == "sector, year, age 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 sector, age, 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 == "sector") %>%   
  filter (interaction == "none") %>%
  mutate (estimate = (exp(estimate) - 1) * 100) %>%  
  select (ssyk, estimate, statistic.y, interaction) %>%
  rename (`F-value` = statistic.y) %>%
  rename (`Increase in salary` = estimate) %>%
  arrange (desc (`F-value`)) %>%
  knitr::kable(
    booktabs = TRUE,
    caption = 'Correlation for F-value (sector) 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 == "sector and sex") %>%
  mutate (estimate = (exp(estimate) - 1) * 100) %>%  
  select (ssyk, estimate, statistic.y, interaction) %>%
  rename (`F-value` = statistic.y) %>%
  rename (`Increase in salary` = estimate) %>%
  arrange (desc (`F-value`)) %>%
  knitr::kable(
    booktabs = TRUE,
    caption = 'Correlation for F-value (sector 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 == "sector and age") %>%
  mutate (estimate = (exp(estimate) - 1) * 100) %>%  
  select (ssyk, estimate, statistic.y, interaction) %>%
  rename (`F-value` = statistic.y) %>%
  rename (`Increase in salary` = estimate) %>%
  arrange (desc (`F-value`)) %>%
  knitr::kable(
    booktabs = TRUE,
    caption = 'Correlation for F-value (sector and age) 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 == "sector and year") %>%
  mutate (estimate = (exp(estimate) - 1) * 100) %>%  
  select (ssyk, estimate, statistic.y, interaction) %>%
  rename (`F-value` = statistic.y) %>%
  rename (`Increase in salary` = estimate) %>%
  arrange (desc (`F-value`)) %>%
  knitr::kable(
    booktabs = TRUE,
    caption = 'Correlation for F-value (sector 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 == "sector, year, age and sex") %>%
  filter (!(contcol.y < 3 & interaction == "sector, year, age and sex")) %>%  
  mutate (estimate = (exp(estimate) - 1) * 100) %>%  
  select (ssyk, estimate, statistic.y, interaction) %>%
  rename (`F-value` = statistic.y) %>%
  rename (`Increase in salary` = estimate) %>%
  arrange (desc (`F-value`)) %>%
  knitr::kable(
    booktabs = TRUE,
    caption = 'Correlation for F-value (sector, year, age and sex) and the yearly increase in salaries')

Let’s check what we have found.

temp <- tb %>%
  filter(`occuptional  (SSYK 2012)` == "218 Specialists within environmental and health protection")
 
model <- lm (log(salary) ~ year_n + poly(age_n, 3) + sector + sex, data = temp) 
 
plot_model(model, type = "pred", terms = c("age_n", "sector"))
## 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 sector, Specialists within environmental and health protection

temp <- tb %>%
  filter(`occuptional  (SSYK 2012)` == "512 Cooks and cold-buffet managers")
 
model <-lm (log(salary) ~ year_n + poly(age_n, 3) + sector + sex, data = temp) 
 
plot_model(model, type = "pred", terms = c("age_n", "sector"))
## Model has log-transformed response. Back-transforming predictions to original response scale. Standard errors are still on the log-scale.

Figure 4: Lowest F-value sector, Cooks and cold-buffet managers

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

Figure 5: Highest F-value interaction sector and gender, Cleaners and helpers

temp <- tb %>%
  filter(`occuptional  (SSYK 2012)` == "511 Cabin crew, guides and related workers")
 
model <- lm (log(salary) ~ year_n + poly(age_n, 3) + sector * sex, data = temp) 
 
plot_model(model, type = "pred", terms = c("age_n", "sex", "sector"))
## Model has log-transformed response. Back-transforming predictions to original response scale. Standard errors are still on the log-scale.

Figure 6: Lowest F-value interaction sector and gender, Cabin crew, guides and related workers

temp <- tb %>%
  filter(`occuptional  (SSYK 2012)` == "251 ICT architects, systems analysts and test managers")
 
model <- lm (log(salary) ~ year_n + poly(age_n, 3) * sector + sex, data = temp) 
 
plot_model(model, type = "pred", terms = c("age_n", "sector"))
## Model has log-transformed response. Back-transforming predictions to original response scale. Standard errors are still on the log-scale.

Figure 7: Highest F-value interaction sector and age, ICT architects, systems analysts and test managers

temp <- tb %>%
  filter(`occuptional  (SSYK 2012)` == "611 Market gardeners and crop growers")
 
model <- lm (log(salary) ~ year_n + poly(age_n, 3) * sector + sex, data = temp) 
 
plot_model(model, type = "pred", terms = c("age_n", "sector"))
## Model has log-transformed response. Back-transforming predictions to original response scale. Standard errors are still on the log-scale.

Figure 8: Lowest F-value interaction sector and age, Market gardeners and crop growers

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

Figure 9: Highest F-value interaction sector and year, Nursing professionals

temp <- tb %>%
  filter(`occuptional  (SSYK 2012)` == "218 Specialists within environmental and health protection")
 
model <- lm (log(salary) ~ year_n * sector + poly(age_n, 3) + sex, data = temp) 
 
plot_model(model, type = "pred", terms = c("age_n", "year_n", "sector"))
## Model has log-transformed response. Back-transforming predictions to original response scale. Standard errors are still on the log-scale.

Figure 10: Lowest F-value interaction sector and year, Specialists within environmental and health protection

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

Figure 11: Highest F-value interaction sector, age, year and gender, Engineering professionals

## TableGrob (2 x 1) "arrange": 2 grobs
##   z     cells    name           grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (2-2,1-1) arrange gtable[layout]
temp <- tb %>%
  filter(`occuptional  (SSYK 2012)` == "243 Marketing and public relations professionals")
 
model <- lm (log(salary) ~ year_n * poly(age_n, 3) * sector * sex, data = temp) 
 
plot_model(model, type = "pred", terms = c("age_n", "year_n", "sex", "sector"))
## 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 sector, age, year and gender, Marketing and public relations professionals

## TableGrob (2 x 1) "arrange": 2 grobs
##   z     cells    name           grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (2-2,1-1) arrange gtable[layout]