Per cent who are women in different occupational groups in Sweden, feature importance

In a previous post, I analysed the feature importance for the per cent of engineers in Sweden who are women. I found that the size of the region is a feature that is significant for the per cent of engineers in Sweden who are women. In this post, I will analyse the feature importance of different occupational groups in Sweden. I will use an ensemble of linear models in my analysis.

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

Please send suggestions for improvement of the analysis to ranalystatisticssweden@gmail.com.

First, define libraries and functions.

library (tidyverse)
## -- Attaching packages -------------------------------------------------- tidyverse 1.3.0 --
## v ggplot2 3.3.0     v purrr   0.3.4
## v tibble  3.0.0     v dplyr   0.8.5
## v tidyr   1.0.2     v stringr 1.4.0
## v readr   1.3.1     v forcats 0.5.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 (caret)    
## Loading required package: lattice
## 
## Attaching package: 'caret'
## The following object is masked from 'package:purrr':
## 
##     lift
library (recipes)  
## 
## Attaching package: 'recipes'
## The following object is masked from 'package:stringr':
## 
##     fixed
## The following object is masked from 'package:stats':
## 
##     step
library (PerformanceAnalytics)
## Loading required package: xts
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## 
## Attaching package: 'xts'
## The following objects are masked from 'package:dplyr':
## 
##     first, last
## 
## Attaching package: 'PerformanceAnalytics'
## The following object is masked from 'package:graphics':
## 
##     legend
library (ggpubr)
## Loading required package: magrittr
## 
## Attaching package: 'magrittr'
## The following object is masked from 'package:purrr':
## 
##     set_names
## The following object is masked from 'package:tidyr':
## 
##     extract
library (ipred) 
library (iml)
library (SuperLearner)
## Loading required package: nnls
## Super Learner
## Version: 2.0-26
## Package created on 2019-10-27
library (scatterplot3d)

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 = groupsize) %>%
  drop_na() %>%
  mutate (year_n = parse_number (year))
}

perc_women <- function(x){  
  ifelse (length(x) == 2, x[2] / (x[1] + x[2]), NA)
} 

nuts <- read.csv("nuts.csv") %>%
  mutate(NUTS2_sh = substr(NUTS2, 3, 4))

nuts %>% 
  distinct (NUTS2_en) %>%
  knitr::kable(
    booktabs = TRUE,
    caption = 'Nomenclature des Unités Territoriales Statistiques (NUTS)')

SL.lm.caret <- function(..., method = "lm", tuneLength = 3, obsWeights = obsWeights, trControl = caret::trainControl(method = "cv", number = 10, verboseIter = FALSE)){
    suppressWarnings(SL.caret(..., obsWeights = obsWeights, method = method, tuneLength = tuneLength, trControl = trControl))
}

SL.lmStepAIC.caret <- function(..., method = "lmStepAIC", tuneLength = 3, obsWeights = obsWeights, trControl = caret::trainControl(method = "cv", number = 10, verboseIter = FALSE)){
    suppressWarnings(SL.caret(..., obsWeights = obsWeights, method = method, tuneLength = tuneLength, trControl = trControl))
}  

SL.bayesglm.caret <- function(..., method = "bayesglm", tuneLength = 3, obsWeights = obsWeights, trControl = caret::trainControl(method = "cv", number = 10, verboseIter = FALSE)){
    suppressWarnings(SL.caret(..., obsWeights = obsWeights, method = method, tuneLength = tuneLength, trControl = trControl))
}  

SL.rlm.caret <- function(..., method = "rlm", tuneLength = 3, obsWeights = obsWeights, trControl = caret::trainControl(method = "cv", number = 10, verboseIter = FALSE)){
    suppressWarnings(SL.caret(..., obsWeights = obsWeights, method = method, tuneLength = tuneLength, trControl = trControl))
}

The data tables are downloaded from Statistics Sweden. They are saved as a comma-delimited file without heading, UF0506A1.csv, http://www.statistikdatabasen.scb.se/pxweb/en/ssd/.

The tables:

UF0506A1_1.csv: Population 16-74 years of age by region, highest level of education, age and sex. Year 1985 – 2018 NUTS 2 level 2008- 10 year intervals (16-74)

000000CG_1: 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.

000000CD_1.csv: 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 Number of employees All sectors.

The data is aggregated, the size of each group is in the column groupsize.

I have also included some calculated predictors from the original data.

perc_women: The percentage of women within each group defined by edulevel, region and year

perc_women_region: The percentage of women within each group defined by year and region

regioneduyears: The average number of education years per capita within each group defined by year and region

eduquotient: The quotient between regioneduyears for men and women

salaryquotient: The quotient between salary for men and women within each group defined by year and region

perc_women_eng_region: The percentage of women who are engineers within each group defined by year and region

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

numedulevel[, 2] <- data.frame(c(8, 9, 10, 12, 13, 15, 22, NA))

tb <- readfile("000000CG_1.csv") 
tb <- readfile("000000CD_1.csv") %>% 
  left_join(tb, by = c("region", "year", "sex", "sector","occuptional  (SSYK 2012)")) 

tb <- readfile("UF0506A1_1.csv") %>%  
  right_join(tb, by = c("region", "year", "sex")) %>%
  right_join(numedulevel, by = c("level of education" = "level.of.education")) %>%
  filter(!is.na(eduyears)) %>%  
  mutate(edulevel = `level of education`) %>%
  group_by(edulevel, region, year, sex, `occuptional  (SSYK 2012)`) %>%
  mutate(groupsize_all_ages = sum(groupsize)) %>%  
  group_by(edulevel, region, year, `occuptional  (SSYK 2012)`) %>% 
  mutate (perc_women = perc_women (groupsize_all_ages[1:2])) %>% 
  mutate (suming = sum(groupsize.x)) %>%
  mutate (salary = (groupsize.y[2] * groupsize.x[2] + groupsize.y[1] * groupsize.x[1])/(groupsize.x[2] + groupsize.x[1])) %>%
  group_by (sex, year, region, `occuptional  (SSYK 2012)`) %>%
  mutate(regioneduyears_sex = sum(groupsize * eduyears) / sum(groupsize)) %>%
  mutate(regiongroupsize = sum(groupsize)) %>% 
  mutate(suming_sex = groupsize.x) %>%
  group_by(region, year, `occuptional  (SSYK 2012)`) %>%
  mutate (sum_pop = sum(groupsize)) %>%
  mutate (regioneduyears = sum(groupsize * eduyears) / sum(groupsize)) %>%
  mutate (perc_women_region = perc_women (regiongroupsize[1:2])) %>% 
  mutate (eduquotient = regioneduyears_sex[2] / regioneduyears_sex[1]) %>% 
  mutate (salary_sex = groupsize.y) %>%
  mutate (salaryquotient = salary_sex[2] / salary_sex[1]) %>%   
  mutate (perc_women_eng_region = perc_women(suming_sex[1:2])) %>%  
  left_join(nuts %>% distinct (NUTS2_en, NUTS2_sh), by = c("region" = "NUTS2_en")) %>%
  drop_na()

summary(tb)
##     region              age            level of education     sex           
##  Length:29050       Length:29050       Length:29050       Length:29050      
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##      year             groupsize          year_n        sector         
##  Length:29050       Min.   :   405   Min.   :2014   Length:29050      
##  Class :character   1st Qu.: 25412   1st Qu.:2015   Class :character  
##  Mode  :character   Median : 61291   Median :2016   Mode  :character  
##                     Mean   : 71345   Mean   :2016                     
##                     3rd Qu.:113524   3rd Qu.:2017                     
##                     Max.   :271889   Max.   :2018                     
##  occuptional  (SSYK 2012)  groupsize.x       year_n.x     groupsize.y    
##  Length:29050             Min.   :  100   Min.   :2014   Min.   : 20200  
##  Class :character         1st Qu.:  490   1st Qu.:2015   1st Qu.: 28900  
##  Mode  :character         Median : 1300   Median :2016   Median : 33900  
##                           Mean   : 3258   Mean   :2016   Mean   : 37066  
##                           3rd Qu.: 3400   3rd Qu.:2017   3rd Qu.: 42100  
##                           Max.   :45000   Max.   :2018   Max.   :133600  
##     year_n.y       eduyears       edulevel         groupsize_all_ages
##  Min.   :2014   Min.   : 8.00   Length:29050       Min.   :   405    
##  1st Qu.:2015   1st Qu.: 9.00   Class :character   1st Qu.: 25412    
##  Median :2016   Median :12.00   Mode  :character   Median : 61291    
##  Mean   :2016   Mean   :12.71                      Mean   : 71345    
##  3rd Qu.:2017   3rd Qu.:15.00                      3rd Qu.:113524    
##  Max.   :2018   Max.   :22.00                      Max.   :271889    
##    perc_women         suming          salary       regioneduyears_sex
##  Min.   :0.3575   Min.   :  240   Min.   : 20661   Min.   :11.18     
##  1st Qu.:0.4343   1st Qu.: 1330   1st Qu.: 29046   1st Qu.:11.63     
##  Median :0.4655   Median : 3100   Median : 34041   Median :11.78     
##  Mean   :0.4775   Mean   : 6515   Mean   : 37105   Mean   :11.83     
##  3rd Qu.:0.5132   3rd Qu.: 7400   3rd Qu.: 42068   3rd Qu.:12.09     
##  Max.   :0.6423   Max.   :60000   Max.   :113976   Max.   :12.55     
##  regiongroupsize    suming_sex       sum_pop        regioneduyears 
##  Min.   :128262   Min.   :  100   Min.   : 262870   Min.   :11.39  
##  1st Qu.:292864   1st Qu.:  490   1st Qu.: 596546   1st Qu.:11.56  
##  Median :528643   Median : 1300   Median :1057419   Median :11.82  
##  Mean   :499413   Mean   : 3258   Mean   : 998826   Mean   :11.83  
##  3rd Qu.:708813   3rd Qu.: 3400   3rd Qu.:1417931   3rd Qu.:11.93  
##  Max.   :827940   Max.   :45000   Max.   :1655215   Max.   :12.41  
##  perc_women_region  eduquotient      salary_sex     salaryquotient  
##  Min.   :0.4831    Min.   :1.019   Min.   : 20200   Min.   :0.6423  
##  1st Qu.:0.4890    1st Qu.:1.027   1st Qu.: 28900   1st Qu.:0.9144  
##  Median :0.4937    Median :1.032   Median : 33900   Median :0.9556  
##  Mean   :0.4931    Mean   :1.033   Mean   : 37066   Mean   :0.9502  
##  3rd Qu.:0.4971    3rd Qu.:1.040   3rd Qu.: 42100   3rd Qu.:0.9941  
##  Max.   :0.5014    Max.   :1.047   Max.   :133600   Max.   :1.3090  
##  perc_women_eng_region   NUTS2_sh        
##  Min.   :0.01659       Length:29050      
##  1st Qu.:0.30876       Class :character  
##  Median :0.56000       Mode  :character  
##  Mean   :0.52565                         
##  3rd Qu.:0.72414                         
##  Max.   :0.94527
tbtemp <- ungroup(tb) %>% dplyr::select(salary, suming, year_n, sum_pop, regioneduyears, perc_women_region, salaryquotient, eduquotient, perc_women_eng_region, `occuptional  (SSYK 2012)`)

tb_unique <- unique(tbtemp)

I will use SuperLearner to train the ensemble consisting of four linear models without interactions. The four models are Linear Regression (lm), Linear Regression with Stepwise Selection (lmStepAIC), Bayesian Generalized Linear Model (bayesglm) and Robust Linear Model (rlm).

summary_table = vector()
cor_table = vector()
sp_table <- vector()
rmse_table <- vector()

for (i in unique(tb_unique$`occuptional  (SSYK 2012)`)){
  temp <- filter(tb_unique, `occuptional  (SSYK 2012)` == i)
  if (dim(temp)[1] > 20){
     temp_weights = temp$suming
     temp <- dplyr::select(temp, - c(`occuptional  (SSYK 2012)`, suming))
     blueprint <- recipe(perc_women_eng_region ~ ., data = temp) %>%
       step_integer(matches("Qual|Cond|QC|Qu")) %>%
       step_center(all_numeric(), -all_outcomes()) %>%
       step_scale(all_numeric(), -all_outcomes()) %>%
       step_dummy(all_nominal(), -all_outcomes(), one_hot = TRUE)
     prepare <- prep(blueprint, training = temp)
     temp <- bake(prepare, new_data = temp)
  
     invisible(capture.output(model <- SuperLearner(
        temp$perc_women_eng_region,
        data.frame(dplyr::select(temp, -c(perc_women_eng_region))),
        family = gaussian(),
        verbose = FALSE,
        obsWeights = temp_weights,
        SL.library = list("SL.lm.caret", "SL.lmStepAIC.caret", "SL.bayesglm.caret", "SL.rlm.caret"))))

     pred <- function(object, newdata){
       predict(model, newdata=newdata, onlySL = TRUE)$pred
     }  
    
     predictor <- Predictor$new(model, 
        data = dplyr::select(temp, -perc_women_eng_region), 
        y = temp$perc_women_eng_region,
        predict.fun = pred)   
   
     imp <- FeatureImp$new(predictor, loss = "mae", n.repetitions = 30)
    
     summary_table <- rbind(summary_table, mutate(tibble(.rows = 7), importance = imp$results$importance, feature = imp$results$feature, importance.05 = imp$results$importance.05, ssyk = i))
    
     cor_table <- rbind(cor_table, mutate(tibble(.rows = 7), feature = colnames(dplyr::select(temp, -c(perc_women_eng_region))), cor = cor(dplyr::select(temp, -c(perc_women_eng_region)), temp$perc_women_eng_region), ssyk = i))
    
     sp_table <- rbind(sp_table, mutate(tibble(.rows = 4), coef = model$coef, model = names(model$coef),  ssyk = i))
    
     prs <- postResample(pred = predict(model)$pred, obs = temp$perc_women_eng_region)
    
     rmse_table <- rbind(rmse_table, mutate(tibble(.rows = 1), RMSE = prs[1], Rsquared = prs[2], MAE = prs[3], ssyk = i))    
  }
}
## Registered S3 methods overwritten by 'lme4':
##   method                          from
##   cooks.distance.influence.merMod car 
##   influence.merMod                car 
##   dfbeta.influence.merMod         car 
##   dfbetas.influence.merMod        car

The table below shows the feature values for the different occupation groups and if there is a single important feature (diff1) or if there are two important features (diff2) for the occupational group. The Rsquared value shows if the model for the occupational group does have a good fit.

summary_table %>% 
  group_by(ssyk) %>% 
  group_by(ssyk) %>% 
  dplyr::mutate(diff1 = importance.05[1] / importance[2]) %>% 
  dplyr::mutate(diff2 = importance.05[2] / importance[3]) %>% 
  left_join(cor_table, by = c("ssyk", "feature")) %>% 
  left_join(sp_table %>% spread(model, coef), by=c("ssyk")) %>% 
  left_join(rmse_table, by=c("ssyk")) %>% 
  dplyr::select(ssyk, feature, importance, importance.05, diff1, diff2, Rsquared) %>%
  knitr::kable( 
     booktabs = TRUE,
     caption = 'Feature values for different occupation groups')

The sum of the per cent that the model was used by the SuperLearner analysing the different occupational groups.

sp_table %>%
  ggplot (aes(coef, model)) +  
    geom_col ()  

Figure 1: The sum of the per cent that the model was used by the SuperLearner

The sum of the strongest feature for every occupational group.

summary_table %>% 
  arrange(desc(importance)) %>% 
  group_by(ssyk) %>% 
  slice(1) %>%
  ggplot (aes(importance, feature)) +  
    geom_col () 

Figure 2: The sum of the strongest feature for every occupational group

Let’s see what we have found. First, check the occupation groups with a single feature that is significantly stronger than all other features. Linear models will not be suitable for all occupational groups implying that the model will not have a high R squared value.

A strong signal, the average number of education years in the region, Personal care workers in health services

temp <- filter(tb_unique, `occuptional  (SSYK 2012)` == "532 Personal care workers in health services")

model <- lm(perc_women_eng_region ~ regioneduyears, weights = suming, data = temp)

temp %>%
  ggplot () +  
    geom_jitter (mapping = aes(x = regioneduyears, y = perc_women_eng_region, colour = suming)) +
    geom_abline (slope = model$coefficients[2], intercept = model$coefficients[1])  +
    labs(
      x = "Education years",
      y = "Per cent of women in the occupation"
    )

Figure 3: Personal care workers in health services, Year 2014 – 2018

summary(model)$adj.r.squared
## [1] 0.7732263
anova(model)
## Analysis of Variance Table
## 
## Response: perc_women_eng_region
##                Df  Sum Sq Mean Sq F value    Pr(>F)    
## regioneduyears  1 315.573 315.573  133.98 5.039e-14 ***
## Residuals      38  89.506   2.355                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
postResample(pred = predict(model), obs = temp$perc_women_eng_region)
##       RMSE   Rsquared        MAE 
## 0.01225219 0.69069055 0.01023249

A strong signal, the average number of education years in the region, Medical doctors

temp <- filter(tb_unique, `occuptional  (SSYK 2012)` == "221 Medical doctors")

model <- lm(perc_women_eng_region ~ regioneduyears, weights = suming, data = temp)

temp %>%
  ggplot () +  
    geom_jitter (mapping = aes(x = regioneduyears, y = perc_women_eng_region, colour = suming)) +
    geom_abline(slope = model$coefficients[2], intercept = model$coefficients[1]) +
    labs(
      x = "Education years",
      y = "Per cent of women in the occupation"
    )

Figure 4: Medical doctors, Year 2014 – 2018

summary(model)$adj.r.squared
## [1] 0.8057127
anova(model)
## Analysis of Variance Table
## 
## Response: perc_women_eng_region
##                Df  Sum Sq Mean Sq F value    Pr(>F)    
## regioneduyears  1 164.765 164.765  154.44 1.385e-14 ***
## Residuals      36  38.407   1.067                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
postResample(pred = predict(model), obs = temp$perc_women_eng_region)
##       RMSE   Rsquared        MAE 
## 0.01683530 0.72088034 0.01385548

A strong signal, the per cent women in the region, Insurance advisers, sales and purchasing agents

temp <- filter(tb_unique, `occuptional  (SSYK 2012)` == "332 Insurance advisers, sales and purchasing agents")

model <- lm(perc_women_eng_region ~ perc_women_region, weights = suming, data = temp)

temp %>%
  ggplot () +  
    geom_jitter (mapping = aes(x = perc_women_region, y = perc_women_eng_region, colour = suming)) +
    geom_abline(slope = model$coefficients[2], intercept = model$coefficients[1]) +
    labs(
      x = "Per cent of women in the region",
      y = "Per cent of women in the occupation"
    )

Figure 5: Insurance advisers, sales and purchasing agents, Year 2014 – 2018

summary(model)$adj.r.squared
## [1] 0.6283407
anova(model)
## Analysis of Variance Table
## 
## Response: perc_women_eng_region
##                   Df Sum Sq Mean Sq F value    Pr(>F)    
## perc_women_region  1 529.66  529.66  56.791 1.395e-08 ***
## Residuals         32 298.45    9.33                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
postResample(pred = predict(model), obs = temp$perc_women_eng_region)
##       RMSE   Rsquared        MAE 
## 0.02935038 0.49206133 0.02250770

Two strong signals, population size in the region and the average number of education years in the region, Engineering professionals

temp <- filter(tb_unique, `occuptional  (SSYK 2012)` == "214 Engineering professionals")

s3d <- scatterplot3d(
  temp$sum_pop, 
  temp$regioneduyears, 
  temp$perc_women_eng_region,
  type = "h", 
  color = "blue", 
  xlab = "Population in region",
  ylab = "Education years",
  zlab = "Per cent of women in the occupation")

model <- lm(perc_women_eng_region ~ sum_pop + regioneduyears, weights = suming, data = temp)

s3d$plane3d(model)

Figure 6: Engineering professionals, Year 2014 – 2018

summary(model)$adj.r.squared
## [1] 0.8121964
anova(model)
## Analysis of Variance Table
## 
## Response: perc_women_eng_region
##                Df  Sum Sq Mean Sq F value    Pr(>F)    
## sum_pop         1 255.902 255.902 144.321 5.673e-14 ***
## regioneduyears  1  31.373  31.373  17.693 0.0001712 ***
## Residuals      35  62.060   1.773                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
postResample(pred = predict(model), obs = temp$perc_women_eng_region)
##        RMSE    Rsquared         MAE 
## 0.012229213 0.835386966 0.009935413

Two strong signals, population size in the region and the per cent women in the region, Insurance advisers, sales and purchasing agents

temp <- filter(tb_unique, `occuptional  (SSYK 2012)` == "332 Insurance advisers, sales and purchasing agents")

s3d <- scatterplot3d(
  temp$sum_pop, 
  temp$perc_women_region, 
  temp$perc_women_eng_region,
  type = "h", 
  color = "blue", 
  xlab = "Population in region",
  ylab = "Per cent of women in the region",
  zlab = "Per cent of women in the occupation")

model <- lm(perc_women_eng_region ~ sum_pop + perc_women_region, weights = suming, data = temp)

s3d$plane3d(model)

Figure 7: Insurance advisers, sales and purchasing agents, Year 2014 – 2018

summary(model)$adj.r.squared
## [1] 0.6525952
anova(model)
## Analysis of Variance Table
## 
## Response: perc_women_eng_region
##                   Df Sum Sq Mean Sq F value    Pr(>F)    
## sum_pop            1 263.40 263.403  30.214 5.168e-06 ***
## perc_women_region  1 294.45 294.455  33.776 2.099e-06 ***
## Residuals         31 270.25   8.718                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
postResample(pred = predict(model), obs = temp$perc_women_eng_region)
##       RMSE   Rsquared        MAE 
## 0.02638844 0.57325855 0.02034915

Two strong signals, year and the per cent women in the region, Physical and engineering science technicians

temp <- filter(tb_unique, `occuptional  (SSYK 2012)` == "311 Physical and engineering science technicians")

s3d <- scatterplot3d(
  temp$year_n, 
  temp$perc_women_region, 
  temp$perc_women_eng_region, 
  type = "h", 
  color = "blue", 
  xlab = "Year",
  ylab = "Per cent of women in the region",
  zlab = "Per cent of women in the occupation")

model <- lm(perc_women_eng_region ~ year_n + perc_women_region, weights = suming, data = temp)

s3d$plane3d(model)

Figure 8: Physical and engineering science technicians, Year 2014 – 2018

summary(model)$adj.r.squared
## [1] 0.5373011
anova(model)
## Analysis of Variance Table
## 
## Response: perc_women_eng_region
##                   Df Sum Sq Mean Sq F value    Pr(>F)    
## year_n             1  32.63  32.630  7.6503  0.009621 ** 
## perc_women_region  1 134.39 134.393 31.5091 4.127e-06 ***
## Residuals         30 127.96   4.265                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
postResample(pred = predict(model), obs = temp$perc_women_eng_region)
##       RMSE   Rsquared        MAE 
## 0.01695193 0.59082239 0.01266243

Two strong signals, year and salary, Naprapaths, physiotherapists, occupational therapists

temp <- filter(tb_unique, `occuptional  (SSYK 2012)` == "227 Naprapaths, physiotherapists, occupational therapists")

s3d <- scatterplot3d(
  temp$year_n, 
  temp$salary, 
  temp$perc_women_eng_region, 
  type = "h", 
  color = "blue", 
  xlab = "Year",
  ylab = "Salary",
  zlab = "Per cent of women in the occupation")

model <- lm(perc_women_eng_region ~ year_n + salary, weights = suming, data = temp)

s3d$plane3d(model)

Figure 9: Naprapaths, physiotherapists, occupational therapists, Year 2014 – 2018

summary(model)$adj.r.squared
## [1] 0.5269917
anova(model)
## Analysis of Variance Table
## 
## Response: perc_women_eng_region
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## year_n     1 5.8240  5.8240  16.077 0.0005492 ***
## salary     1 4.9902  4.9902  13.776 0.0011481 ** 
## Residuals 23 8.3317  0.3622                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
postResample(pred = predict(model), obs = temp$perc_women_eng_region)
##       RMSE   Rsquared        MAE 
## 0.01261698 0.46523146 0.01003402