The significance of education on the salary of engineers in Sweden

In my last posts, I analysed the significance of experience for different occupational groups. In this post, I will turn the interest towards education. I will again start with engineers and see if I can expand my analysis to all occupational groups.

First, define libraries and functions.

library (tidyverse) 
## -- Attaching packages -------------------------------------------- tidyverse 1.2.1 --
## v ggplot2 3.2.0     v purrr   0.3.2
## v tibble  2.1.3     v dplyr   0.8.3
## v tidyr   0.8.3     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 (splines)
#install_github("ZheyuanLi/SplinesUtils")
library(SplinesUtils)

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, 000000CY.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 sector, occupational group (SSYK 2012), sex and educational level (SUN). Year 2014 – 2018 Monthly salary All sectors

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

The column level of education is renamed because TukeyHSD doesn’t handle variable names within quotes.

tb <- readfile("000000CY.csv") %>% 
  filter(`occuptional  (SSYK 2012)` == "214 Engineering professionals") %>%       
  mutate(edulevel = `level of education`)

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

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

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

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

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

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

Figure 2: 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)$adj.r.squared  
## [1] 0.8327314
Anova(model, type = 2) %>% 
  tidy() %>% 
  knitr::kable( 
  booktabs = TRUE,
  caption = 'Anova report from linear model fit')

How much do the different levels of education affect the salary? We can calculate the differences between the levels with Tukey’s honest significant difference. All significant level differences are shown in the table below.

tukeytable <- TukeyHSD(aov(log(salary) ~ sex + edulevel, data = tb)) %>%
  tidy() %>%
  mutate(percdiff = (1 / exp(estimate) - 1) * 100)

tukeytable  %>% 
  filter(adj.p.value < 0.05) %>%
  arrange(estimate) %>%  
  knitr::kable( 
  booktabs = TRUE,
  caption = 'Tukey HSD 95 % confidence intervals for the pairwise significant differences')

We can conclude from the summary table that there is a positive correlation between longer education and higher salaries.

From the table of Tukey’s honest significant difference, we can see the difference in salaries between the different education lengths. Note that the estimates are negative due to the log transformation, the untransformed differences are in the column percdiff.

Can we approximate how much the salaries increase by one year of education by assigning a numeric value to the factors in the table?

As a first approach, I will use the data in the following table.

I will use a B-spline function to approximate the increase in salaries over age. The rows for “no information about the level of educational attainment” is removed from the table from Statistics Sweden.

numedulevel <- read.csv("edulevel.csv") 

numedulevel %>%
  knitr::kable(
  booktabs = TRUE,
  caption = 'Initial approach, length of education') 

tbnum <- tb %>% 
    right_join(numedulevel, by = c("level of education" = "level.of.education")) %>%
  filter(!is.na(eduyears))
## Warning: Column `level of education`/`level.of.education` joining character
## vector and factor, coercing into character vector
modelcont <- lm(log(salary) ~ bs(eduyears, knots = c(14)) + year_n + sex, data = tbnum)

tbnum <- bind_cols(tbnum, as_tibble(exp(predict(modelcont, tbnum, interval = "confidence"))))

The F-value from the Anova table for years is 146 and the adjusted R-squared value is 0,932 implying a good fit of the model. Both the F-value and the adjusted R-squared increased from then using the categorical predictors. (Removing the rows with “no information about level of educational attainment” improves the adjusted R-squared for the model with categorical predictors to 0.931.)

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

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

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

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

What does the continous function from the model look like?

contspline <- RegBsplineAsPiecePoly(modelcont, "bs(eduyears, knots = c(14))")

tibble(eduyears = 9:19) %>%
  ggplot () + 
    geom_point (mapping = aes(x = eduyears,y = predict(contspline, eduyears))) +
  labs(
    x = "Years of education",
    y = "Salary"
  )

Figure 4: Model fit, SSYK 214, Correlation between education and salary

And it’s derivative.

tibble(eduyears = 9:19) %>%
  ggplot () + 
    geom_point (mapping = aes(x = eduyears,y = (exp(predict(contspline, eduyears, deriv = 1)) - 1) * 100)) +
  labs(
    x = "Years of education",
    y = "Salary difference (%)"
  )

Figure 5: Model fit, SSYK 214, The derivative for education

Comparison between the categorical and the continuous predictor. Column withinconf states if the estimate from the numerical model is within the 95 % confidence interval from the Tukey’s honest significant difference table. All estimates from the model with the continuous predictor are within the 95 % confidence intervals from the Tukey HSD table.

tukeytable <- tukeytable %>% 
  rowwise() %>% 
  mutate(comp_from = unlist(strsplit(comparison, ")-"))[1]) %>%
  rowwise() %>% 
  mutate(comp_to = unlist(strsplit(comparison, ")-"))[2]) %>%             mutate(comp_from = paste (comp_from, ")", sep="")) %>%
  left_join(numedulevel, by = c("comp_from" = "level.of.education")) %>%
  left_join(numedulevel, by = c("comp_to" = "level.of.education")) %>%    mutate(numestimate = predict(contspline, eduyears.x) -           predict(contspline, eduyears.y)) %>% 
  mutate(withinconf = numestimate > conf.low && numestimate < conf.high) %>% 
  mutate(percdiffcont = (1 / exp(predict(contspline, eduyears.x) - predict(contspline, eduyears.y)) - 1) * 100)
## Warning: Column `comp_from`/`level.of.education` joining character vector and
## factor, coercing into character vector
## Warning: Column `comp_to`/`level.of.education` joining character vector and
## factor, coercing into character vector
tukeytable %>% 
  select(term, comparison, estimate, adj.p.value, numestimate, withinconf, percdiff, percdiffcont) %>%
  filter(adj.p.value < 0.05) %>%
  arrange(estimate) %>%
  knitr::kable( 
  booktabs = TRUE,
  caption = 'Comparison between categorical and continous predictor')

Now. let´s perform some diagnostics on the models. First, a look at the residuals for the model shows no apparent problem. We can see that the outlier at row 55 has disappeared in the plots for the continuous model.

Three out of four outliers when using categorical predictors were from the factor “no information about level of educational attainment”.

For the continuous predictors there are two outliers in the factor “upper secondary education, 2 years or less (ISCED97 3C)” and two outliers in the factor “upper secondary education 3 years (ISCED97 3A)” indicating that the model could be improved.

tb[20,]$edulevel
## [1] "no information about level of educational attainment"
tb[41,]$edulevel
## [1] "no information about level of educational attainment"
tb[55,]$edulevel
## [1] "no information about level of educational attainment"
tb[57,]$edulevel
## [1] "upper secondary education, 2 years or less (ISCED97 3C)"
tbnum[12,]$edulevel
## [1] "upper secondary education, 2 years or less (ISCED97 3C)"
tbnum[18,]$edulevel
## [1] "upper secondary education, 2 years or less (ISCED97 3C)"
tbnum[25,]$edulevel
## [1] "upper secondary education 3 years (ISCED97 3A)"
tbnum[29,]$edulevel
## [1] "upper secondary education 3 years (ISCED97 3A)"
plot(model, which = 1)

Figure 6: Residuals vs Fitted of model fit

plot(modelcont, which = 1)

Figure 7: Residuals vs Fitted of model fit

The Normal Q-Q shows some possible outliers.

plot(model, which = 2)

Figure 8: Normal Q-Q

plot(modelcont, which = 2)

Figure 9: Normal Q-Q

Again, the Standardised residuals show some possible outliers.

plot(model, which = 3)

Figure 10: Scale-Location

plot(modelcont, which = 3)

Figure 11: Scale-Location

The outliers are also found in the Leverage plot.

plot(model, which = 5)

Figure 12: Residuals vs Leverage

plot(modelcont, which = 5)

Figure 13: Residuals vs Leverage