In my previous posts, I have shown that there exists a correlation between salaries and experience as well as between salaries and education within many occupational groups. Within the data I used from Statistics Sweden there was also information about gender. Some of the models and graphs also showed that there is a difference in salaries between the different genders.

In this post, I will examine the interaction between the predictors. Apart from gender the interaction between year also needs to be investigated.

It would have been very interesting to examine the interaction between experience and education. However, this information is not available at Statistics Sweden, at least not for free.

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   1.0.0     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)
library(sjPlot)

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, 00000031.csv, http://www.statistikdatabasen.scb.se/pxweb/en/ssd/.

Average monthly pay (total pay), non-manual workers private sector (SLP), SEK by occupational group (SSYK), age, sex and year. Year 2014 – 2018

There are no women in the age group 65-66 years which explains the uncertainty for this group.

We are interested to see if the interaction terms between gender, age and year are important factors in salaries. As a null hypothesis, we assume that the interaction terms between gender, age and year is not related to the salary and examine if we can reject this hypothesis with the data from Statistics Sweden.

tb <- readfile("00000031.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) 
   
tb %>%
  ggplot () +  
    geom_point (mapping = aes(x = year_n,y = salary, colour = age, 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 + age, data = tb)

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

summary(model)$adj.r.squared  
## [1] 0.9817476
Anova(model, type = 2)
## Anova Table (Type II tests)
## 
## Response: log(salary)
##            Sum Sq Df F value    Pr(>F)    
## year_n    0.06878  1 117.504 < 2.2e-16 ***
## sex       0.02581  1  44.095 3.166e-09 ***
## age       2.81366  9 534.074 < 2.2e-16 ***
## Residuals 0.04800 82                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(model, which = 1)

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

tb[38,]
## Source: local data frame [1 x 9]
## Groups: <by row>
## 
## # A tibble: 1 x 9
##   `occuptional  (SSYK 2012~ age      sex   year  salary year_n age_l age_h age_n
##   <chr>                     <chr>    <chr> <chr>  <dbl>  <dbl> <int> <int> <dbl>
## 1 214 Engineering professi~ 18-24 y~ men   2016   27300   2016    18    24    21
tb[56,]
## Source: local data frame [1 x 9]
## Groups: <by row>
## 
## # A tibble: 1 x 9
##   `occuptional  (SSYK 2012~ age      sex   year  salary year_n age_l age_h age_n
##   <chr>                     <chr>    <chr> <chr>  <dbl>  <dbl> <int> <int> <dbl>
## 1 214 Engineering professi~ 65-66 y~ men   2016   51500   2016    65    66  65.5
tb[72,]
## Source: local data frame [1 x 9]
## Groups: <by row>
## 
## # A tibble: 1 x 9
##   `occuptional  (SSYK 2012~ age      sex   year  salary year_n age_l age_h age_n
##   <chr>                     <chr>    <chr> <chr>  <dbl>  <dbl> <int> <int> <dbl>
## 1 214 Engineering professi~ 55-59 y~ women 2017   45000   2017    55    59    57

I will examine models with both the categorical and continuous predictors to see if they can make the same predictions. The models with continuous predictors show approximately the same as the models with categorical predictors as you can see below.

We start by examining the interaction between year and age. We find that the only significant interaction between age and year is in the age group 65-66 years old.

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

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

summary(model)$adj.r.squared  
## [1] 0.9817432
Anova(model, type = 2)
## Anova Table (Type II tests)
## 
## Response: log(salary)
##             Sum Sq Df  F value    Pr(>F)    
## sex        0.02581  1  44.0847 4.817e-09 ***
## year_n     0.06878  1 117.4762 < 2.2e-16 ***
## age        2.81366  9 533.9464 < 2.2e-16 ***
## year_n:age 0.00526  9   0.9978    0.4497    
## Residuals  0.04274 73                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot_model(model, type = "pred", terms = c("year_n", "age"), colors = "Set3")
## Model has log-transformed response. Back-transforming predictions to original response scale. Standard errors are still on the log-scale.

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

plot(model, which = 1)

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

model <- lm (log(salary) ~ sex + year_n * bs(age_n, knots = c(35, 50)), data = tb)

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

summary(model)$adj.r.squared  
## [1] 0.98315
Anova(model, type = 2)
## Anova Table (Type II tests)
## 
## Response: log(salary)
##                                      Sum Sq Df   F value    Pr(>F)    
## sex                                 0.02618  1   48.4439 8.049e-10 ***
## year_n                              0.06883  1  127.3747 < 2.2e-16 ***
## bs(age_n, knots = c(35, 50))        2.81344  5 1041.2680 < 2.2e-16 ***
## year_n:bs(age_n, knots = c(35, 50)) 0.00445  5    1.6454    0.1574    
## Residuals                           0.04377 81                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot_model(model, type = "pred", terms = c("age_n", "year_n"))
## 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

Examine the interaction between gender and age. We can see that women have a significantly lower salary for the age groups 35-39 years, 45-49 years, 50-54 years, 55-59 years and 60-64 years.

The F-value from the Anova table for the interaction term between gender and age is 3,2 (Pr(>F) < 0.003564), sufficient for rejecting the null hypothesis that the interaction between gender and age has no effect on the salary holding year as constant.

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

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

summary(model)$adj.r.squared  
## [1] 0.9849736
Anova(model, type = 2)
## Note: model has aliased coefficients
##       sums of squares computed by model comparison
## Anova Table (Type II tests)
## 
## Response: log(salary)
##            Sum Sq Df  F value    Pr(>F)    
## year_n    0.06878  1 142.7313 < 2.2e-16 ***
## sex       0.02581  1  53.5620 2.492e-10 ***
## age       2.81366  9 648.7344 < 2.2e-16 ***
## sex:age   0.01234  8   3.2006  0.003564 ** 
## Residuals 0.03566 74                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot_model(model, type = "pred", terms = c("age", "sex"))
## Warning in predict.lm(model, newdata = fitfram, type = "response", se.fit =
## se, : prediction from a rank-deficient fit may be misleading
## 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[75,]
## Source: local data frame [1 x 9]
## Groups: <by row>
## 
## # A tibble: 1 x 9
##   `occuptional  (SSYK 2012~ age      sex   year  salary year_n age_l age_h age_n
##   <chr>                     <chr>    <chr> <chr>  <dbl>  <dbl> <int> <int> <dbl>
## 1 214 Engineering professi~ 65-66 y~ men   2017   46500   2017    65    66  65.5
model <- lm (log(salary) ~ year_n + sex * bs(age_n, knots = c(35, 50)), data = tb)

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

summary(model)$adj.r.squared  
## [1] 0.9859063
Anova(model, type = 2)
## Anova Table (Type II tests)
## 
## Response: log(salary)
##                                   Sum Sq Df   F value    Pr(>F)    
## year_n                           0.06880  1  152.2042 < 2.2e-16 ***
## sex                              0.02612  1   57.7814 4.558e-11 ***
## bs(age_n, knots = c(35, 50))     2.81344  5 1244.9030 < 2.2e-16 ***
## sex:bs(age_n, knots = c(35, 50)) 0.01161  5    5.1354 0.0003844 ***
## Residuals                        0.03661 81                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot_model(model, type = "pred", terms = c("age_n", "sex"))
## 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

Examine the interaction between gender and year. Both the model with categorical predictors and the model with continuous predictors indicates that women’s salaries increase more than men’s salaries.

The F-value from the Anova table for the interaction term between gender and year is 3,1 (Pr(>F) < 0.08053), with 90 % significance we can reject the null hypothesis that the interaction between gender and year has no effect on the salary holding age as constant.

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

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

summary(model)$adj.r.squared  
## [1] 0.9822101
Anova(model, type = 2)
## Anova Table (Type II tests)
## 
## Response: log(salary)
##             Sum Sq Df F value    Pr(>F)    
## year_n     0.06878  1 120.559 < 2.2e-16 ***
## sex        0.02581  1  45.242 2.266e-09 ***
## age        2.81430  9 548.084 < 2.2e-16 ***
## year_n:sex 0.00179  1   3.132   0.08053 .  
## Residuals  0.04621 81                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot_model(model, type = "pred", terms = c("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

model <- lm (log(salary) ~ year_n * sex + bs(age_n, knots = c(35, 50)), data = tb)

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

summary(model)$adj.r.squared  
## [1] 0.9829655
Anova(model, type = 2)
## Anova Table (Type II tests)
## 
## Response: log(salary)
##                               Sum Sq Df   F value    Pr(>F)    
## year_n                       0.06883  1  125.9950 < 2.2e-16 ***
## sex                          0.02612  1   47.8063 8.182e-10 ***
## bs(age_n, knots = c(35, 50)) 2.81407  5 1030.2198 < 2.2e-16 ***
## year_n:sex                   0.00178  1    3.2607    0.0745 .  
## Residuals                    0.04644 85                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot_model(model, type = "pred", terms = c("year_n", "sex"))
## Model has log-transformed response. Back-transforming predictions to original response scale. Standard errors are still on the log-scale.

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

Examine the interaction between year, gender and age. It’s becoming a bit difficult to make any significant predictions about the interactions from the model, the data is not sufficiently large, the degrees of freedom has been reduced from 82 in the model with no interaction terms to 56.

The following interaction terms are significant to 90 %.

year_n:sexwomen

sexwomen:age35-39 years

year_n:sexwomen:age35-39 years

In the graph below you can that the salaries for men in the age group 35-39 years has increased more than for women.

The F-value from the Anova table for the interaction term between gender and age is 3,3 (Pr(>F) < 0.00373), sufficient for rejecting the null hypothesis that the interaction between gender and age has no effect on the salary holding year as constant.

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

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

summary(model)$adj.r.squared  
## [1] 0.9854223
Anova(model, type = 2)
## Note: model has aliased coefficients
##       sums of squares computed by model comparison
## Anova Table (Type II tests)
## 
## Response: log(salary)
##                 Sum Sq Df  F value    Pr(>F)    
## year_n         0.06878  1 147.1250 < 2.2e-16 ***
## sex            0.02581  1  55.2108 6.727e-10 ***
## age            2.81430  9 668.8557 < 2.2e-16 ***
## year_n:sex     0.00113  1   2.4097   0.12622    
## year_n:age     0.00460  9   1.0927   0.38295    
## sex:age        0.01234  8   3.2991   0.00373 ** 
## year_n:sex:age 0.00310  8   0.8277   0.58192    
## Residuals      0.02618 56                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot_model(model, type = "pred", terms = c("year_n", "age", "sex"))
## Warning in predict.lm(model, newdata = fitfram, type = "response", se.fit =
## se, : prediction from a rank-deficient fit may be misleading
## Model has log-transformed response. Back-transforming predictions to original response scale. Standard errors are still on the log-scale.
## Warning in RColorBrewer::brewer.pal(n, pal): n too large, allowed maximum for palette Set1 is 9
## Returning the palette you asked for with that many colors
## Warning: Removed 10 rows containing missing values (geom_path).

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

plot(model, which = 1)

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

model <- lm (log(salary) ~ year_n * sex * bs(age_n, knots = c(35, 50)), data = tb)

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

summary(model)$adj.r.squared  
## [1] 0.9873689
Anova(model, type = 2)
## Anova Table (Type II tests)
## 
## Response: log(salary)
##                                          Sum Sq Df  F value    Pr(>F)    
## year_n                                  0.06880  1 169.8286 < 2.2e-16 ***
## sex                                     0.03780  6  15.5511 3.011e-11 ***
## bs(age_n, knots = c(35, 50))            2.82569 10 697.5506 < 2.2e-16 ***
## year_n:sex                              0.00126  1   3.1090 0.0822248 .  
## year_n:bs(age_n, knots = c(35, 50))     0.00394  5   1.9452 0.0977419 .  
## sex:bs(age_n, knots = c(35, 50))        0.01163  5   5.7405 0.0001702 ***
## year_n:sex:bs(age_n, knots = c(35, 50)) 0.00253  5   1.2489 0.2958798    
## Residuals                               0.02836 70                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot_model(model, type = "pred", terms = c("age_n", "year_n", "sex"))
## Model has log-transformed response. Back-transforming predictions to original response scale. Standard errors are still on the log-scale.

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

plot(model, which = 1)

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

Average basic salary, monthly salary and women´s salary as a percentage of men´s salary by sector, occupational group (SSYK 2012), gender and educational level (SUN). Year 2014 - 2018 Monthly salary All sectors

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

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

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

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
numedulevel %>%
  knitr::kable(
  booktabs = TRUE,
  caption = 'Initial approach, length 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 17: SSYK 214, Architects, engineers and related professionals, Year 2014 - 2018

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

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

summary(model)$adj.r.squared  
## [1] 0.9311952
Anova(model, type = 2)
## Anova Table (Type II tests)
## 
## Response: log(salary)
##            Sum Sq Df F value    Pr(>F)    
## year_n    0.05936  1  100.27 1.232e-13 ***
## sex       0.08124  1  137.22 4.476e-16 ***
## edulevel  0.34338  5  116.00 < 2.2e-16 ***
## Residuals 0.03019 51                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(model, which = 1)

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

tb[20,]
## # A tibble: 1 x 8
##   sector  `occuptional  (~ sex   `level of educat~ year  salary year_n edulevel 
##   <chr>   <chr>            <chr> <chr>             <chr>  <dbl>  <dbl> <chr>    
## 1 0 all ~ 214 Engineering~ men   no information a~ 2015   48400   2015 no infor~
tb[41,]
## # A tibble: 1 x 8
##   sector  `occuptional  (~ sex   `level of educat~ year  salary year_n edulevel 
##   <chr>   <chr>            <chr> <chr>             <chr>  <dbl>  <dbl> <chr>    
## 1 0 all ~ 214 Engineering~ women no information a~ 2016   39400   2016 no infor~
tb[55,]
## # A tibble: 1 x 8
##   sector  `occuptional  (~ sex   `level of educat~ year  salary year_n edulevel 
##   <chr>   <chr>            <chr> <chr>             <chr>  <dbl>  <dbl> <chr>    
## 1 0 all ~ 214 Engineering~ women no information a~ 2017   35300   2017 no infor~

We start by examining the interaction between year and education. The salaries for the education upper secondary education, 2 years or less increase more than other educations during the time period 2014-2018, however, this can not be shown with 95% significance from this data.

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

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

summary(model)$adj.r.squared  
## [1] 0.9312476
Anova(model, type = 2)
## Anova Table (Type II tests)
## 
## Response: log(salary)
##                  Sum Sq Df  F value    Pr(>F)    
## sex             0.08094  1 136.8227 2.212e-15 ***
## year_n          0.05936  1 100.3448 3.853e-13 ***
## edulevel        0.34338  5 116.0852 < 2.2e-16 ***
## year_n:edulevel 0.00298  5   1.0078     0.424    
## Residuals       0.02721 46                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot_model(model, type = "pred", terms = c("year_n", "edulevel"))
## Model has log-transformed response. Back-transforming predictions to original response scale. Standard errors are still on the log-scale.

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

plot(model, which = 1)

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

tb[10,]
## # A tibble: 1 x 8
##   sector  `occuptional  (~ sex   `level of educat~ year  salary year_n edulevel 
##   <chr>   <chr>            <chr> <chr>             <chr>  <dbl>  <dbl> <chr>    
## 1 0 all ~ 214 Engineering~ women post-secondary e~ 2014   38800   2014 post-sec~
tb[25,]
## # A tibble: 1 x 8
##   sector  `occuptional  (~ sex   `level of educat~ year  salary year_n edulevel 
##   <chr>   <chr>            <chr> <chr>             <chr>  <dbl>  <dbl> <chr>    
## 1 0 all ~ 214 Engineering~ women post-secondary e~ 2015   41200   2015 post-sec~
tb[29,]
## # A tibble: 1 x 8
##   sector  `occuptional  (~ sex   `level of educat~ year  salary year_n edulevel 
##   <chr>   <chr>            <chr> <chr>             <chr>  <dbl>  <dbl> <chr>    
## 1 0 all ~ 214 Engineering~ men   upper secondary ~ 2016   42100   2016 upper se~
model <- lm(log(salary) ~ sex + year_n * bs(eduyears, knots = c(14)), data = tbnum)

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

summary(model)$adj.r.squared  
## [1] 0.9302821
Anova(model, type = 2)
## Anova Table (Type II tests)
## 
## Response: log(salary)
##                                     Sum Sq Df  F value    Pr(>F)    
## sex                                0.08110  1 135.1856 1.453e-15 ***
## year_n                             0.05941  1  99.0272 2.972e-13 ***
## bs(eduyears, knots = c(14))        0.34299  4 142.9343 < 2.2e-16 ***
## year_n:bs(eduyears, knots = c(14)) 0.00179  4   0.7457    0.5657    
## Residuals                          0.02880 48                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot_model(model, type = "pred", terms = c("year_n", "eduyears"))
## Model has log-transformed response. Back-transforming predictions to original response scale. Standard errors are still on the log-scale.

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

plot(model, which = 1)

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

tb[25,]
## # A tibble: 1 x 8
##   sector  `occuptional  (~ sex   `level of educat~ year  salary year_n edulevel 
##   <chr>   <chr>            <chr> <chr>             <chr>  <dbl>  <dbl> <chr>    
## 1 0 all ~ 214 Engineering~ women post-secondary e~ 2015   41200   2015 post-sec~

Examine the interaction between gender and education. The only significant difference in salaries for the length in education is for the group upper secondary education, 2 years or less.

The F-value from the Anova table for the interaction term between gender and education is 3,3 (Pr(>F) < 0.01195), sufficient for rejecting the null hypothesis that the interaction between gender and education has no effect on the salary holding year as constant.

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

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

summary(model)$adj.r.squared  
## [1] 0.9439866
Anova(model, type = 2)
## Anova Table (Type II tests)
## 
## Response: log(salary)
##               Sum Sq Df  F value    Pr(>F)    
## sex          0.08124  1 168.5530 < 2.2e-16 ***
## edulevel     0.34338  5 142.4863 < 2.2e-16 ***
## year_n       0.05879  1 121.9676 1.576e-14 ***
## sex:edulevel 0.00802  5   3.3293   0.01195 *  
## Residuals    0.02217 46                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot_model(model, type = "pred", terms = c("edulevel", "sex"))
## Model has log-transformed response. Back-transforming predictions to original response scale. Standard errors are still on the log-scale.

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

plot(model, which = 1)

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

tb[12,]
## # A tibble: 1 x 8
##   sector  `occuptional  (S~ sex   `level of educa~ year  salary year_n edulevel 
##   <chr>   <chr>             <chr> <chr>            <chr>  <dbl>  <dbl> <chr>    
## 1 0 all ~ 214 Engineering ~ women post-graduate e~ 2014   45000   2014 post-gra~
model <- lm(log(salary) ~ sex * bs(eduyears, knots = c(14)) + year_n, data = tbnum)

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

summary(model)$adj.r.squared  
## [1] 0.943092
Anova(model, type = 2)
## Anova Table (Type II tests)
## 
## Response: log(salary)
##                                  Sum Sq Df  F value    Pr(>F)    
## sex                             0.08128  1 165.9766 < 2.2e-16 ***
## bs(eduyears, knots = c(14))     0.34299  4 175.1086 < 2.2e-16 ***
## year_n                          0.05891  1 120.2975 1.128e-14 ***
## sex:bs(eduyears, knots = c(14)) 0.00708  4   3.6147   0.01181 *  
## Residuals                       0.02350 48                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot_model(model, type = "pred", terms = c("eduyears", "sex"))
## Model has log-transformed response. Back-transforming predictions to original response scale. Standard errors are still on the log-scale.

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

plot(model, which = 1)

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

Examine the interaction between gender and year There is no evidence that the salaries for men and women are interacting with year, holding education constant.

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

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

summary(model)$adj.r.squared  
## [1] 0.93001
Anova(model, type = 2)
## Anova Table (Type II tests)
## 
## Response: log(salary)
##             Sum Sq Df  F value    Pr(>F)    
## sex        0.08124  1 134.8940 8.258e-16 ***
## year_n     0.05936  1  98.5705 2.047e-13 ***
## edulevel   0.34343  5 114.0493 < 2.2e-16 ***
## sex:year_n 0.00008  1   0.1364    0.7134    
## Residuals  0.03011 50                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot_model(model, type = "pred", terms = c("year_n", "sex"))
## Model has log-transformed response. Back-transforming predictions to original response scale. Standard errors are still on the log-scale.

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

plot(model, which = 1)

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

tb[18,]
## # A tibble: 1 x 8
##   sector  `occuptional  (~ sex   `level of educat~ year  salary year_n edulevel 
##   <chr>   <chr>            <chr> <chr>             <chr>  <dbl>  <dbl> <chr>    
## 1 0 all ~ 214 Engineering~ men   post-secondary e~ 2015   44800   2015 post-sec~
model <- lm(log(salary) ~ sex * year_n + bs(eduyears, knots = c(14)), data = tbnum)

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

summary(model)$adj.r.squared  
## [1] 0.9304967
Anova(model, type = 2)
## Anova Table (Type II tests)
## 
## Response: log(salary)
##                              Sum Sq Df  F value    Pr(>F)    
## sex                         0.08128  1 135.8986 5.361e-16 ***
## year_n                      0.05941  1  99.3330 1.446e-13 ***
## bs(eduyears, knots = c(14)) 0.34304  4 143.3973 < 2.2e-16 ***
## sex:year_n                  0.00008  1   0.1401    0.7098    
## Residuals                   0.03050 51                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot_model(model, type = "pred", terms = c("eduyears", "sex"))
## Model has log-transformed response. Back-transforming predictions to original response scale. Standard errors are still on the log-scale.

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

plot(model, which = 1)

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

Examine the interaction between year, gender and education It’s becoming a bit difficult to make any significant predictions from the model, the data is not sufficiently large, the degrees of freedom has been reduced from 51 in the model with no interaction terms to 35.

There are no significant interaction terms in the summary.

The F-value from the Anova table for the interaction term between gender and education is 3,0 (Pr(>F) < 0.02377), sufficient for rejecting the null hypothesis that the interaction between gender and education has no effect on the salary holding year as constant.

model <- lm(log(salary) ~ sex * year_n * edulevel, data = tbnum)

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

summary(model)$adj.r.squared  
## [1] 0.9377879
Anova(model, type = 2)
## Anova Table (Type II tests)
## 
## Response: log(salary)
##                      Sum Sq Df  F value    Pr(>F)    
## sex                 0.08897  6  27.6985 6.428e-12 ***
## year_n              0.05879  1 109.8150 2.456e-12 ***
## edulevel            0.35144 10  65.6498 < 2.2e-16 ***
## sex:year_n          0.00007  1   0.1338   0.71677    
## sex:edulevel        0.00800  5   2.9902   0.02377 *  
## year_n:edulevel     0.00298  5   1.1127   0.37150    
## sex:year_n:edulevel 0.00039  5   0.1445   0.98038    
## Residuals           0.01874 35                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot_model(model, type = "pred", terms = c("year_n", "edulevel","sex"))
## Model has log-transformed response. Back-transforming predictions to original response scale. Standard errors are still on the log-scale.

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

plot(model, which = 1)

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

tb[29,]
## # A tibble: 1 x 8
##   sector  `occuptional  (~ sex   `level of educat~ year  salary year_n edulevel 
##   <chr>   <chr>            <chr> <chr>             <chr>  <dbl>  <dbl> <chr>    
## 1 0 all ~ 214 Engineering~ men   upper secondary ~ 2016   42100   2016 upper se~
model <- lm(log(salary) ~ sex * year_n * bs(eduyears, knots = c(14)), data = tbnum)

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

summary(model)$adj.r.squared  
## [1] 0.9364598
Anova(model, type = 2)
## Anova Table (Type II tests)
## 
## Response: log(salary)
##                                         Sum Sq Df  F value    Pr(>F)    
## sex                                    0.08817  5  32.2521 7.295e-13 ***
## year_n                                 0.05891  1 107.7412 8.809e-13 ***
## bs(eduyears, knots = c(14))            0.35011  8  80.0434 < 2.2e-16 ***
## sex:year_n                             0.00008  1   0.1517   0.69908    
## sex:bs(eduyears, knots = c(14))        0.00706  4   3.2268   0.02218 *  
## year_n:bs(eduyears, knots = c(14))     0.00179  4   0.8168   0.52231    
## sex:year_n:bs(eduyears, knots = c(14)) 0.00032  4   0.1468   0.96335    
## Residuals                              0.02132 39                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot_model(model, type = "pred", terms = c("eduyears", "year_n","sex"))
## Model has log-transformed response. Back-transforming predictions to original response scale. Standard errors are still on the log-scale.

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

plot(model, which = 1)

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