# What color is your paycheck?

Posted on Fri 04 December 2015 in blog

One of the data science skills I want to play around with is deriving insights from data that publically available. Here, lets use some data on SF employee compensation and see what we can learn from the data.

First, per usual, load the dependencies.

```
%matplotlib inline
import requests
import json
import pandas as pd
import matplotlib.pylab as plt
import numpy as np
import requests
import seaborn as sns
import sys
sys.path.append('../../dataSandbox/forPelican/') #for loading custom plotting script
from sfBokeh import plotScatter, plotPCA
```

### Get the data from the API in json format¶

```
URL = "https://data.sfgov.org/resource/88g8-5mnd.json"
response = requests.get(URL)
response.raise_for_status()
data = response.json()
```

### Put it into a pandas dataframe¶

Notice that for some downstream steps, we should convert the data to float. Then, take a look at the data

```
df = pd.DataFrame.from_dict(data)
quantVar = ['health_dental','other_benefits','other_salaries','overtime','retirement','salaries','total_benefits','total_compensation','total_salary']
df[quantVar] = df[quantVar].astype(float)
df = df.sort_values(by='total_compensation')
df = df.reset_index()
df.head()
```

### A first look¶

There are a lot of quantitative and qualitative variables (columns) per employee (row). This should be a great setup for an unsupervised learning approach. But I don't want to give anything away yet...

Look at one variable and you can see some relationships between the variables, shown in the equations below.

```
df.loc[600]
```

### Exploring salary and benefit combinations¶

Let's make a bar plot sorted for total employee compensation to see how total compenation is broken down in terms of salary and benefits

```
fig = plt.figure(figsize=(5,4), dpi=1600)
ax = plt.subplot(111)
annoIdx = set(['total_salary','total_benefits'])
colors = ['red','blue']
width = 1
condn = range(0,len(df.index))
rollBottom = np.repeat(0,len(condn)) #Create an index so the bars are stacked on top of each other
condnIdx = np.arange(len(condn)) # the x locations for the groups
colIdx = 0
for i in annoIdx:
ax.bar(condnIdx, df[i], width, color=colors[colIdx],bottom=rollBottom,alpha=0.2,edgecolor = "none")
rollBottom = rollBottom + list(df[i])
colIdx = colIdx + 1
box = ax.get_position()
ax.set_position([box.x0, box.y0, box.width * 0.8, box.height*0.95])
# Put a legend to the right of the current axis
ax.legend([i for i in annoIdx],loc='center left', bbox_to_anchor=(1, 0.5))
ax.set_xlabel('Sorted employee index',fontsize=16,fontweight='bold')
ax.set_ylabel('Total compensation',fontsize=16,fontweight='bold')
plt.show()
```

### Show me the money¶

Unsurprisingly, salary is the big contributor to total compensation. There are a number of departments in the data. Who makes the most? What is the spread? Let's use pandas with a seaborn boxplot to take a look.

```
ob = df.groupby(['department','department_code']).mean()[['other_salaries','total_compensation','overtime']]
tcOrder = [i for i in ob.sort_values(by='total_compensation',ascending=False).index.get_level_values(1)]
fig = plt.figure(figsize=(10,4), dpi=1600)
ax = sns.boxplot(x = 'department_code',y='total_compensation',data=df,order=tcOrder)#,color='skyblue')#,order=otOrder)
labs = [i.get_text() for i in ax.get_xticklabels()]
ax.set_xticklabels(labs,rotation=45)
ax.set_xlabel('Department code',fontsize=16)
ax.set_ylabel('Total compensation',fontsize=16)
plt.show()
```

### Digging deeper¶

So city attorneys (CAT) have a sweet gig apparently. Let's look at the individual variables underlying the total salary/benefits. Pandas has some great functionality here, and making the data points transparent reveals data density.

```
fig = plt.figure(figsize=(10,9), dpi=1600)
q2 = ['health_dental','other_benefits','other_salaries','overtime','retirement','salaries','total_benefits','total_compensation','total_salary','department']
ax = pd.tools.plotting.scatter_matrix(df[quantVar],figsize=(10, 10),alpha=0.2)
#sns.pairplot(df[q2],hue='department')
plt.show()
```

### Follow the variance¶

There are a lot of things going on here. Rather than trying to understand each plot one-by-one, let's use Principal Component Analysis (PCA). PCA has several benefits. It identifies

- the axes of greatest variance in the dataset (which may be more important)
- lower dimensional reconstructions of the data to ease interpretability and relationships between variables.

We will use the `sklearn`

package to do the PCA. Also, we will remove `total_compensation`

, `total_benefits`

, and `total_salary`

, since we explicitly know these are related and PCA works best with independent variables. Also, let's scale the data, so that the total salary doesn't just wash out the signal.

```
from sklearn.preprocessing import scale
from sklearn.decomposition import PCA
fig = plt.figure(figsize=(6,4), dpi=1600)
ax = plt.subplot(111)
dfs = pd.DataFrame(scale(df[quantVar],with_std=True),columns=df[quantVar].columns)
dfResp = dfs['total_compensation']
#Remove as much collinearilty as possible
pred = [quantVar[i] for i in range(len(quantVar)) if quantVar[i] not in ['total_compensation','total_salary','total_benefits']]
dfPred = dfs[pred]
pca = PCA()
pca.fit(dfPred)
expVar = pca.explained_variance_ratio_
pca.explained_variance_ratio_
sns.barplot(x=np.arange(1,len(expVar)+1),y=expVar*100,alpha=0.8,color='darkblue')
ax.set_ylabel('Perc. var. explained',fontsize=16,fontweight='bold')
ax.set_xlabel('Principal Component',fontsize=16,fontweight='bold')
plt.show()
```

### Transform the data¶

Great, it looks like the first 3 PCs explain most of the variance. So now, let's express the data in terms of our shiny new PCs.

```
#New data in PC space
newData = pd.DataFrame(pca.transform(dfPred),columns=['PC' + str(i+1) for i in range(len(expVar))])
#Combine with original data to keep annotations, etc
newData = pd.concat([newData, df], axis=1)
```

### Colorfully and interactively separated¶

Let's first do a plot between PC1 and PC2, the two axes of greatest variance. Also, we can encode data annotations by color to learn what drives differences between departments.

Hover over the points to see the department, job title, and a couple quantitative variables. I chose parks and rec (REC) and the fire department (FIR).

```
plotPCA(myData=newData,dpt1='REC', dpt2='FIR', xVar='PC1', yVar='PC2',legLoc='bottom_left')
```