Salaries in the city of Norfolk municipality

A case study of hierarchical Bayesian modeling.

xxxxxxxxxx
 
md"# Salaries in the city of Norfolk municipality
​
_A case study of hierarchical Bayesian modeling._
"

1.7 ms

xxxxxxxxxx
 
begin
    using CSV, DataFrames, Turing, StatsPlots
    gr()
end;

21.1 s

Data

The data is a list of employees. We consider here the salary as the random variable of interest and department and employee status as explanatory variables.

xxxxxxxxxx
 
md"## Data
​
The data is a list of employees. We consider here the salary as the random variable of interest and department and employee status as explanatory variables.
"

12.7 μs

	Department	Position Title	Employee Status	Initial Hire Date	Date in Position	Fair Labor Standards Act (FLSA)	Base Salary
	String	String	String	String	String	String	Float64
1	"CF - MacArthur Memorial"	"Museum Attendant"	"Casual Part-time"	"06/04/2007"	"06/04/2007"	"Nonexempt"	12.51
2	"CF - MacArthur Memorial"	"Museum Attendant"	"Casual Part-time"	"11/13/2006"	"11/13/2006"	"Nonexempt"	12.51
3	"CF - MacArthur Memorial"	"Museum Attendant"	"Permanent Full-time"	"09/23/2020"	"09/23/2020"	"Nonexempt"	26200.0
4	"CF - MacArthur Memorial"	"Museum Attendant"	"Permanent Full-time"	"01/07/2019"	"01/07/2019"	"Nonexempt"	26200.0
5	"CF - MacArthur Memorial"	"Administrative Technician"	"Permanent Full-time"	"09/23/2015"	"03/19/2018"	"Nonexempt"	33454.1
6	"CF - MacArthur Memorial"	"Curator"	"Permanent Full-time"	"11/03/2014"	"11/03/2014"	"Exempt"	55193.6
7	"CF - MacArthur Memorial"	"Education Manager"	"Permanent Full-time"	"10/29/2007"	"11/23/2009"	"Exempt"	59534.6
8	"CF - MacArthur Memorial"	"Archivist"	"Permanent Full-time"	"06/01/1994"	"05/15/1996"	"Exempt"	61718.1
9	"CF-Cultural&Convention Center"	"Ticket Seller"	"Intermittent Temporary"	"12/30/2004"	"12/30/2004"	"Nonexempt"	11.91
10	"CF-Cultural&Convention Center"	"Ticket Sales Supervisor"	"Intermittent Temporary"	"10/18/2002"	"10/18/2002"	"Nonexempt"	14.59
more
4399	"Zoo-Veterinary & Wellness Camp"	"Veterinarian"	"Permanent Full-time"	"09/26/2020"	"09/26/2020"	"Exempt"	75000.0

xxxxxxxxxx
 
df = CSV.File("02norfolk_employee_data.csv") |> DataFrame

17.6 s

Some of the salaries are per hour, while others are per year. Fortunately, the highest hourly rate is much lower than the lowest yearly salary, hence we can ‘fix’ the data by scaling up hourly rates by the number of work hours per year (≈ 2000).

We find the threshold separating hourly and yearly salaries as the geometric mean of the two ends of the longest interval between subsequent sorted salaries:

xxxxxxxxxx
 
    md"Some of the salaries are per hour, while others are per year. Fortunately, the highest hourly rate is much lower than the lowest yearly salary, hence we can ‘fix’ the data by scaling up hourly rates by the number of work hours per year (≈ 2000).
​
We find the threshold separating hourly and yearly salaries as the geometric mean of the two ends of the longest interval between subsequent sorted salaries:
"

6.6 μs

731.0526656814815

xxxxxxxxxx
 
begin
    salary = df[!, Symbol("Base Salary")]
    department = df[!, Symbol("Department  ")]
    status = df[!, Symbol("Employee Status")]
    threshold = let salary = copy(salary)
        sort!(salary)
        d = salary[2:end] ./ salary[begin:end-1]
        imax = findmax(d)[2]
        threshold = sqrt(salary[imax]*salary[imax+1])
    end 
end

185 ms

xxxxxxxxxx
 
let
    histogram(log.(salary[salary .< threshold]), type=:density, label="hourly", color="lightblue", xlabel="log salary")
    histogram!(log.(salary[salary .>= threshold]), type=:density, label="yearly", color="lightgreen")
    vline!(log.([threshold]), label="threshold", lw=2, color="red")
end

26.0 s

With hourly salaries multiplied by the number of work hours in a year, the distribution of salaries is unimodel.

xxxxxxxxxx
 
md"With hourly salaries multiplied by the number of work hours in a year, the distribution of salaries is unimodel."

6.1 μs

xxxxxxxxxx
 
yearly_salary = let 
    hours_per_year = 40*52 # 40 hours per week, 52 weeks
    yearly = copy(salary)
    yearly[yearly .< threshold] .*= hours_per_year
    yearly
end;

170 ms

xxxxxxxxxx
 
histogram(log.(yearly_salary), normalize=:pdf, xlabel="log salary", label="yearly salary")

210 ms

One can see from the plots that the distribution of all salaries looks more like a Gauss' bell on the log scale. Because of that, we log-transform the data. Afterwards, we standardize the data to zero mean and unit variance.

xxxxxxxxxx
 
md"One can see from the plots that the distribution of all salaries looks more like a Gauss' bell on the log scale. Because of that, we log-transform the data. Afterwards, we standardize the data to zero mean and unit variance.
"

3.4 μs

xxxxxxxxxx
 
begin
    log_yearly_salary = log.(yearly_salary);
    log_yearly_salary_mean = mean(log_yearly_salary)
    log_yearly_salary_std = std(log_yearly_salary)
    log_yearly_salary = (log_yearly_salary .- log_yearly_salary_mean) ./
                        log_yearly_salary_std
end;

49.2 ms

For convenience, the data is transformed into three arrays:

salary –- the target variable,
department –- an explanatory variable,
status –- employee status, another explanatory variable.

We collect the data in two forms: with dependence only on the department, and on both the department and the employee status.

xxxxxxxxxx
 
md"For convenience, the data is transformed into three arrays:
​
* salary --- the target variable,
* department --- an explanatory variable,
* status --- employee status, another explanatory variable.
​
We collect the data in two forms: with dependence only on the department, and on both the department and the employee status.
"

2.8 ms

xxxxxxxxxx
 
begin
    function list2map(list)
        map = Dict{String, Int}()
        for (i, n) in enumerate(unique(list))
            map[n] = i
        end
        map
    end
    deptmap = list2map(department)
    statmap = list2map(status)
end;

38.9 ms

xxxxxxxxxx
 
begin
    deptidx = map((d) -> deptmap[d], department)
    statidx = map((s) -> statmap[s], status)
end;

84.3 μs

Enter cell code...xxxxxxxxxx
 
​

97.0 ns

Models

We define two models, model_by_dept and model_by_dept_stat. model_by_dept is conditioned only on the department. model_by_dept_stat is conditioned on both the department and on employee status.

xxxxxxxxxx
 
md"## Models
​
We define two models, `model_by_dept` and `model_by_dept_stat`. `model_by_dept` is conditioned only on the department. `model_by_dept_stat` is conditioned on both the department and on employee status.
"

30.6 μs

Model by department only

To speed up the inference, instead of conditioning on each employee individually, we collect mean, variance, and count for each department. For a normal distribution, empirical mean $m$ and empirical variance $s^{2}$ are distributed as (see Wikipedia):

$m \sim N o r m a l (μ, \frac{σ}{\sqrt{n}})$

$s^{2} \sim G a m m a (\frac{n - 1}{2}, \frac{2 σ^{2}}{n - 1})$

Conditioning on individual employees would work too but take significantly more time — we would have $\approx 2000$ observation points instead of just $165$ .

xxxxxxxxxx
 
md"### Model by department only
​
To speed up the inference, instead of conditioning on each employee individually, we collect mean, variance, and count for each department. For a normal distribution, empirical mean $m$ and empirical variance $s^2$ are distributed as (see [Wikipedia](https://en.wikipedia.org/wiki/Normal_distribution)):
​
$$m \sim \mathrm{Normal}\left(\mu, \frac \sigma {\sqrt n}\right)$$
$$s^2 \sim \mathrm{Gamma}\left(\frac {n-1} 2, \frac {2\sigma^2} {n-1}\right)$$
​
​
Conditioning on individual employees would work too but take significantly more time — we would have $\approx 2000$ observation points instead of just $165$.
"

18.6 μs

xxxxxxxxxx
 
begin
    by_dept_mean = Array{Float64, 1}(undef, length(deptmap))
    by_dept_var = Array{Float64, 1}(undef, length(deptmap))
    by_dept_count = Array{Float64, 1}(undef, length(deptmap))
    for i in 1:length(deptmap)
        by_dept = log_yearly_salary[deptidx .== i]
        by_dept_mean[i] = mean(by_dept)
        by_dept_var[i] = var(by_dept)
        by_dept_count[i] = length(by_dept)
    end
end

1.4 ms

Conditioning on variance can cause problems of two kinds. First, for less than 2 samples, the variance is undefined. Then, if all salaries are the same, the variance is zero, and the density of Gamma distribution on 0 for k > 1 is 0.

Therefore, we only condition on variance if there are enough samples (more than 3) and the values of the samples are different. This constraint is an indirect consequence of our use of LogNormal approximation of the distribution of salaries.

xxxxxxxxxx
 
md"Conditioning on variance can cause problems of two kinds. First, for less than 2 samples, the variance is undefined. Then, if all salaries are the same, the variance is zero, and the density of Gamma distribution on 0 for k > 1 is 0. 
​
​
Therefore, we only condition on variance if there are enough samples (more than 3) and the values of the samples are different. This constraint is an indirect consequence of our use of LogNormal approximation of the distribution of salaries.
"

8.5 μs

model_by_dept (generic function with 1 method)

xxxxxxxxxx
 
@model function model_by_dept(mean, var, count)
    μμ ~ Normal(0, 1)
    σμ ~ Gamma(1, 1)
    μσ ~ Normal(0, 1)
    σσ ~ Gamma(1, 1)
    μ ~ MvNormal(fill(μμ, length(mean)), σμ)
    σ ~ MvLogNormal(fill(μσ, length(mean)), σσ)
    for i in eachindex(mean)
        σ2 = σ[i]^2
        mean[i] ~ Normal(μ[i], sqrt(σ2/count[i]))
        k = (count[i]-1)/2 
        if k > 1 && var[i] > 0 
            var[i] ~ Gamma(k, σ2/k)
        end
    end
end

204 μs