Lecture 22: Principal Components Analysis – Data 100, Fall 2020¶
Notebook originally by Josh Hug (Fall 2019)
Edits by Anthony D. Joseph and Suraj Rampure (Fall 2020)
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
import numpy as np
We often create visualizations in order to facilitate exploratory data analysis. For example, we might create scatterplots to explore the relationship between pairs of variables in a dataset.
The dataset below gives the "percentage body fat, age, weight, height, and ten body circumference measurements" for 252 men.
http://jse.amstat.org/v4n1/datasets.johnson.html
For simplicity, we read in only 8 of the provided attributes, yielding the given dataframe.
#http://jse.amstat.org/datasets/fat.txt
df3 = pd.read_fwf("data/fat.dat.txt", colspecs = [(17, 21), (23, 29), (35, 37),
(39, 45), (48, 53), (73, 77),
(80, 85), (88, 93)], header=None, names = ["% fat", "density", "age", "weight", "height", "neck", "chest", "abdomen"])
df3.head()
We see that percentage fat and density in g/cm^3 are almost completely redundant.
sns.scatterplot(data = df3, x = "% fat", y = "density");
By contrast, while there is a strong correlation between neck and chest measurements, the resulting data is very 2 dimensional.
sns.scatterplot(data = df3, x = "neck", y = "chest");
Age and height show a small correlation as peolpe seem to get slightly smaller with greater age in this dataset.
sns.scatterplot(data = df3, x = "age", y = "height");
We note that there is one outlier where a person is slightly less than 29.5 inches tall. While there are a extraordinarily small number of adult males who are less than three feet tall, reflection on the rest of the data from this observation suggest that this was simply an error.
df3.query("height < 40")
df3 = df3.drop(41)
sns.scatterplot(data = df3, x = "age", y = "height");
We can try to visualize more than 2 attributes at once, but the relationships displayed in e.g. the color and dot size space are much harder for human readers to see. For example, above we saw that density and % fat are almost entirely redundant, but this relationship is impossible to see when comparing the colors and dot sizes.
sns.scatterplot(data = df3, x = "neck", y = "chest", hue="density", size = "% fat");
Seaborn gives us the ability to create a matrix of all possible pairs of variables. This is can be useful, though even with only 8 variables it's still difficult to fully digest.
sns.pairplot(df3);
We should note that despite the very strong relationship between % fat and density, the numerical rank of the data matrix is still 8. For the rank to be 7, we'd need the data to be almost exactly on a line. We'll talk about techniques to reduce the dimensionality over the course of this lecture and the next.
np.linalg.matrix_rank(df3)
House of Representatives Voting Data¶
Next, let's consider voting data from the house of representatives in the U.S. during the month of September 2019. In this example, our goal will be to try to find clusters of representatives who vote in similar ways. For example, we might expect to find that Democrats and Republicans vote similarly to other members of their party.
from pathlib import Path
from ds100_utils import fetch_and_cache
from datetime import datetime
from IPython.display import display
import yaml
plt.rcParams['figure.figsize'] = (4, 4)
plt.rcParams['figure.dpi'] = 150
sns.set()
base_url = 'https://github.com/unitedstates/congress-legislators/raw/master/'
legislators_path = 'legislators-current.yaml'
f = fetch_and_cache(base_url + legislators_path, legislators_path)
legislators_data = yaml.load(open(f))
def to_date(s):
return datetime.strptime(s, '%Y-%m-%d')
legs = pd.DataFrame(
columns=['leg_id', 'first', 'last', 'gender', 'state', 'chamber', 'party', 'birthday'],
data=[[x['id']['bioguide'],
x['name']['first'],
x['name']['last'],
x['bio']['gender'],
x['terms'][-1]['state'],
x['terms'][-1]['type'],
x['terms'][-1]['party'],
to_date(x['bio']['birthday'])] for x in legislators_data])
legs.head(3)
# February 2019 House of Representatives roll call votes
# Downloaded using https://github.com/eyeseast/propublica-congress
votes = pd.read_csv('data/votes.csv')
votes = votes.astype({"roll call": str})
votes.head()
votes.merge(legs, left_on='member', right_on='leg_id').sample(5)
def was_yes(s):
if s.iloc[0] == 'Yes':
return 1
else:
return 0
vote_pivot = votes.pivot_table(index='member',
columns='roll call',
values='vote',
aggfunc=was_yes,
fill_value=0)
print(vote_pivot.shape)
vote_pivot.head()
vote_pivot.shape
This data has 441 observations (members of the House of Representatives including the 6 non-voting representatives) and 41 dimensions (votes). While politics is quite polarized, none of these columns are linearly dependent as we note below.
np.linalg.matrix_rank(vote_pivot)
Suppose we want to find clusters of similar voting behavior. We might try by reducing our data to only two dimensions and looking to see if we can identify clear patterns. Let's start by looking at what votes were most controversial.
np.var(vote_pivot, axis=0).sort_values(ascending = False)
We see that roll call 548 had very little variance. According to http://clerk.house.gov/evs/2019/roll548.xml, this bill was referring to the 2019 Whistleblower Complaint about President Trump and Ukraine. The full text of the house resolution for this roll call can be found at https://www.congress.gov/bill/116th-congress/house-resolution/576/text:
(1) the whistleblower complaint received on August 12, 2019, by the Inspector General of the Intelligence Community shall be transmitted immediately to the Select Committee on Intelligence of the Senate and the Permanent Select Committee on Intelligence of the House of Representatives; and
(2) the Select Committee on Intelligence of the Senate and the Permanent Select Committee on Intelligence of the House of Representatives should be allowed to evaluate the complaint in a deliberate and bipartisan manner consistent with applicable statutes and processes in order to safeguard classified and sensitive information.
We see that 421 congresspeople voted for this resolution, and 12 did not vote for this resolution. 2 members answered "present" but did not vote no, and 10 did not vote at all. Clearly, a scatterplot involving this particular dimension isn't going to be useful.
vote_pivot['548'].value_counts()
By contrast, we saw high variance for most of the other roll call votes. Most them had variances near 0.25, which is the maximum possible for a variable which can take on values 0 or 1. Let's consider the two highest variance variables, shown below:
vote_pivot['555'].value_counts()
vote_pivot['530'].value_counts()
Let's use these as our two dimensions for our scatterplot and see what happens.
sns.scatterplot(x='530', y='555', data=vote_pivot);
By adding some random noise, we can get rid of the overplotting.
vote_pivot_jittered = vote_pivot.copy()
vote_pivot_jittered.loc[:, '515':'555'] += np.random.random(vote_pivot_jittered.loc[:, '515':'555'].shape) * 0.3
sns.scatterplot(x='530', y='555', data=vote_pivot_jittered);
We can also look at this data labeled by party.
vote_pivot_labeled = vote_pivot.reset_index().merge(legs, left_on='member', right_on='leg_id').set_index('member')
vote_pivot_labeled.head(5)
vote_pivot_labeled_jittered = vote_pivot_labeled.copy()
vote_pivot_labeled_jittered.loc[:, '515':'555'] += np.random.random(vote_pivot_labeled_jittered.loc[:, '515':'555'].shape) * 0.3
sns.scatterplot(x='530', y='555', data=vote_pivot_labeled_jittered, hue="party");
sns.scatterplot(x='545', y='552', data=vote_pivot_labeled_jittered, hue="party");
We see that considering only two votes does seem to do a pretty good job of telling Republicans from Democrats. We'll see later in this lecture how we can do even better using a technique called "Principle Component Analysis" or PCA.
Before we can get there, we'll need to spend some time reviewing key linear algebra principles.
Matrix Operations in Python¶
age_and_height = np.array([[182, 28], [399, 30], [725, 33]])
M = np.array([[1, 0, 0], [0, 1, 1/12]])
age_and_height @ M
Singular Value Decomposition Experiment¶
Manual Decomposition¶
In the table below, we have the width, height, area, and perimeter of a rectangle stored in a dataframe.
rectangle = pd.read_csv("data/rectangle_data.csv")
rectangle.head(5)
Naturally the perimeter is just the sum of 2x the width and 2x the height. Thus, if we create a new dataframe that has only the width, height, and area...
rectangle_no_perim = rectangle[["width", "height", "area"]]
rectangle_no_perim.head(5)
... then we can recover the perimeter by multiplying this matrix by
1 0 0 2
0 1 0 2
0 0 1 0
transform_3D_to_4D = [[1, 0, 0, 2], [0, 1, 0, 2], [0, 0, 1, 0]]
rectangle_with_perimeter_back = np.array(rectangle_no_perim) @ transform_3D_to_4D
pd.DataFrame(rectangle_with_perimeter_back, columns = ["width", "height", "area", "perimeter"]).head(5)
Singular Value Decomposition Example¶
Singular value decomposition is a numerical technique to (among other things) automatically uncover such redundancies. Given an input matrix X, SVD will return $U\Sigma$ and $V^T$ such that $ X = U \Sigma V^T $.
u, s, vt = np.linalg.svd(rectangle, full_matrices = False)
As we did before with our manual decomposition, we can recover our original rectangle data by multiplying the three return values of this function back together.
pd.DataFrame(u * s @ vt).head(4)
The two key pieces of the decomposition are $U\Sigma$ and $V^T$, which we can think of for now as analogous to our 'data' and 'transformation operation' from our manual decomposition earlier.
Let's start by looking at $U\Sigma$, which we can compute with the Python code u*s.
u*s
Similarly, we can look at vt.
vt
The automatic decomposition returned by the svd function looks quite different than what we got when we manually decomposed our data into "data" and "operations". That is, vt is a bunch of seemingly arbitrary numbers instead of the rather simple:
1 0 0 2
0 1 0 2
0 0 1 0
Similarly, if we look at the shape of $U\Sigma$ and $V^T$ we see that they are bigger than in our manual decomposition. Specifically $U\Sigma$ still has 4 columns, meaning that each observation is 4 dimensional. Furthermore, rather than our transformation operation $V^T$ being 3x4, it's 4x4 rows tall, meaning that it maps 4 dimensional inputs back to 4 dimensions.
This seems problematic, because our goal of using SVD was to find a transformation operation that takes 3D inputs and maps them up to 4 dimensions.
Luckily, if we look carefully at $U\Sigma$, we see that the last attribute of each observation is very close to 0.
u * s
Thus, it makes sense that we remove the last column of $U \Sigma$.
u = u[:, 0:3]
s = s[0:3]
u * s
Similarly, because the observations are now 3D, we should remove the last row of $V^T$, since we want to use $V^T$ to map our now 3D data into 4D (rather than expecting 4D input data).
vt = vt[0:3, :]
After removing the redundant portions of $U\Sigma$ and $V^T$, we can verify that multiplying them together again yields our original array.
u * s @ vt
Singular Value Decomposition Experiment¶
The reasons that $U \Sigma$ and $V^T$ look so different than the results of our manual decomposition are a consequence of how singular value decomposition. Specifically, given X, SVD computers $U$, $\Sigma$, and $V$ such that:
- $\Sigma$ is a diagonal matrix containing the singular values of X.
- $U$ and $V$ are matrices whose columns form an orthonormal set.
- $X = U \Sigma V^T$
That is, there are an infinite number of matrices such that $X = AB$. In our example above, we created A and B through manual calculation and insight (recognizing that perimeter was 2w + 2h). SVD computes them automatically, and the results have the specific properties enumerated above.
u, s, vt = np.linalg.svd(rectangle, full_matrices = False)
The Singular Values Matrix $\Sigma$¶
The middle result returned by the SVD process is a diagonal matrix consisting of the singular values. We won't say a lot about what they mean today, but we will note that if the matrix is of rank r, then the first r singular values will be non-zero, and the rest will be zero.
s
Python returns the singular values in a slightly funny format (as a list). To get them into the correct form, we can use diag(s).
np.diag(s)
If we do this, then we can compute the original matrix using only the matrix multiplication operator. That is, rather than writing u * s @ vt, we can do u @ np.diag(s) @ vt.
(u @ np.diag(s) @ vt)[0:5, :]
The fact that the last singular value is zero is why we were able to remove the last dimension of $U$ and the last operation in $V^T$. That is, since the last column of $U \Sigma$ is always zero, we can just remove it completely.
(u @ np.diag(s))
(u[:, :-1] @ np.diag(s[:-1]) @ vt[:-1, :])[0:5, :]
$U$ and $V^T$¶
Let's try verifying that $U$ is orthonormal. If so, any column dot producted with itself should be 1, and any columns dot producted with any other should be zero.
np.dot(u[:, 0], u[:, 0])
np.dot(u[:, 1], u[:, 1])
np.dot(u[:, 2], u[:, 2])
for i in range(0, u.shape[1]):
print(f"row {i} dot producted with itself is {np.dot(u[:, i], u[:, i])}")
for i in range(0, u.shape[1]):
for j in range(i + 1, u.shape[1]):
print(f"row {i} dot producted with {j} is {np.dot(u[:, i], u[:, j])}")
Let's now look at $V$.
v = vt.T
for i in range(0, v.shape[1]):
print(f"row {i} of V dot producted with itself is {np.dot(v[:, i], v[:, i])}")
for i in range(0, v.shape[1]):
for j in range(i + 1, v.shape[1]):
print(f"row {i} of V dot producted with {j} is {np.dot(v[:, i], v[:, j])}")
We can also see that the transpose of u (and v) is also the inverse of u (and v).
u.T @ u
v.T @ v
Low Rank Approximations of Rectangle Data¶
This section was not covered in lecture but may add additional context.
colnames = list(rectangle.columns)
Suppose we our area and perimeter measurements are slightly noisy.
rectangle_noisy = rectangle.copy()
rectangle_noisy["perimeter"] += np.random.normal(loc = 0, scale = 0.01, size = len(rectangle_noisy["perimeter"]))
rectangle_noisy["area"] += np.random.normal(loc = 0, scale = 0.01, size = len(rectangle_noisy["area"]))
rectangle_noisy.head(5)
u, s, vt = np.linalg.svd(rectangle_noisy, full_matrices = False)
Even on this noisy data, we can recover the exact measurements by multiplying $U \Sigma V^T$.
vt
pd.DataFrame((u @ np.diag(s) @ vt)[0:5, :]).head(5)
Let's take a look at the singular values:
s
Despite not being zero, we can nevertheless again eliminate the last column (dimension) of $U$, last singular value from $\Sigma$ and last row (operation) from $V^T$. We see that we get back reasonable answers, though they're a little off.
pd.DataFrame((u[:, :-1] @ np.diag(s[:-1]) @ vt[:-1, :]),
columns = colnames).head(5)
This is what is known as a "low rank approximation" of the data. Specifically, this is a rank 3 approximation since we're using 3 of the available 4 dimensions. Since the noise was small, this rank 3 approximation is good.
As a warmup to the rest of lecture, let's see what the rank 2 approximation looks like.
Note: In the code below, I've specified which rows/columns that I do want rather than the ones I don't want. That is, instead of writing u[:, :-2], I've written u[:, 0:2].
pd.DataFrame((u[:, 0:2] @ np.diag(s[0:2]) @ vt[0:2, :]),
columns = colnames).head(5)
This code is getting a little hard to read, so let's write a function that computes the rank k approximation of a given DataFrame.
def compute_rank_k_approximation(data, k):
u, s, vt = np.linalg.svd(data, full_matrices = False)
return pd.DataFrame(u[:, 0:k] @ np.diag(s[0:k]) @ vt[0:k, :], columns = data.columns)
pd.DataFrame(compute_rank_k_approximation(rectangle_noisy, 2)).head(5)
We can even look at the rank 1 approximation.
compute_rank_k_approximation(rectangle_noisy, 1).head(5)
In effect, what the SVD process has done is figured out a way to rebuild an approximation of our entire 4D matrix using only 1 column of data. This might seem somewhat shocking.
pd.DataFrame((u[:, 0:1] @ np.diag(s[0:1]) @ vt[0:1, :]),
columns = colnames).head(5)
u, s, vt = np.linalg.svd(rectangle, full_matrices = False)
compute_rank_k_approximation(rectangle_noisy, 1).head(5)
def compute_rank_k_approximation(data, k):
u, s, vt = np.linalg.svd(data, full_matrices = False)
return pd.DataFrame(u[:, 0:k] @ np.diag(s[0:k]) @ vt[0:k, :], columns = data.columns)
Rank 1 Approximation of 2D Data, Data Centering¶
# Downloads from https://www.gapminder.org/data/
cm_path = 'data/child_mortality_0_5_year_olds_dying_per_1000_born.csv'
fe_path = 'data/children_per_woman_total_fertility.csv'
cm = pd.read_csv(cm_path).set_index('country')['2017'].to_frame()/10
fe = pd.read_csv(fe_path).set_index('country')['2017'].to_frame()
child_data = cm.merge(fe, left_index=True, right_index=True).dropna()
child_data.columns = ['mortality', 'fertility']
child_data.head()
def scatter14(data):
sns.scatterplot('mortality', 'fertility', data=data)
plt.xlim([0, 14])
plt.ylim([0, 14])
plt.xticks(np.arange(0, 14, 2))
plt.yticks(np.arange(0, 14, 2))
scatter14(child_data)
sns.scatterplot('mortality', 'fertility', data=child_data)
u, s, vt = np.linalg.svd(child_data, full_matrices = False)
child_data_reconstructed = pd.DataFrame(u @ np.diag(s) @ vt, columns = ["mortality", "fertility"], index=child_data.index)
As we'd expect, the product of $U$, $\Sigma$, and $V^T$ recovers the original data perfectly.
child_data_reconstructed.head(5)
sns.scatterplot('mortality', 'fertility', data=child_data)
What happens if we throw away a column of $U$, a singular value from $\Sigma$, and a row from $V^T$? In this case we end up with the "rank 1 approximation" of the data.
Looking at the data, we see that it does a surprisingly good job.
#child_data_rank_1_approximation = pd.DataFrame(u[:, :-1] @ np.diag(s[:-1]) @ vt[:-1, :], columns = ["mortality", "fertility"], index=child_data.index)
child_data_rank_1_approximation = compute_rank_k_approximation(child_data, 1)
child_data_rank_1_approximation.head(5)
By plotting the data in a 2D space, we can see what's going on. We're simply getting the original data projected on to some 1 dimensional subspace.
sns.scatterplot('mortality', 'fertility', data=child_data_rank_1_approximation)
There's one significant issue with our projection, which we can see by plotting both the original data and our reconstruction on the same axis. The issue is that the projection goes through the origin but our data has a non-zero y-intercept.
sns.scatterplot('mortality', 'fertility', data=child_data)
sns.scatterplot('mortality', 'fertility', data=child_data_rank_1_approximation)
While this y-intercept misalignment isn't terrible here, it can be really bad. For example, consider the 2D dataset below (from our body measurements dataset from the previous lecture).
#http://jse.amstat.org/datasets/fat.txt
body_data = pd.read_fwf("data/fat.dat.txt", colspecs = [(9, 13), (17, 21), (23, 29), (35, 37),
(39, 45), (48, 53), (57, 61), (64, 69),
(73, 77), (80, 85), (88, 93), (96, 101),
(105, 109), (113, 117), (121, 125), (129, 133),
(137, 141), (145, 149)],
header=None, names = ["% brozek fat", "% siri fat", "density", "age",
"weight", "height", "adiposity", "fat free weight",
"neck", "chest", "abdomen", "hip", "thigh",
"knee", "ankle", "bicep", "forearm",
"wrist"])
#body_data = body_data.drop(41) #drop the weird record
body_data.head()
density_and_abdomen = body_data[["density", "abdomen"]]
density_and_abdomen.head(5)
If we look at the data, the rank 1 approximation looks at least vaguely sane from the table.
density_and_abdomen_rank_1_approximation = compute_rank_k_approximation(density_and_abdomen, 1)
density_and_abdomen_rank_1_approximation.head(5)
But if we plot on 2D axes, we'll see that things are very wrong.
sns.scatterplot(x="density", y="abdomen", data=body_data)
density_and_abdomen_rank_1_approximation = compute_rank_k_approximation(density_and_abdomen, 1)
sns.scatterplot(x="density", y="abdomen", data=body_data)
sns.scatterplot(x="density", y="abdomen", data=density_and_abdomen_rank_1_approximation);
Since the subspace that we're projecting on to is off and to the right, we end up with a bizarre result where our rank 1 approximation believes that density increases with abdomen size, even though the data shows the opposite.
To fix this issue, we should always start the SVD process by zero-centering our data. That is, for each column, we should subtract the mean of that column.
np.mean(density_and_abdomen, axis = 0)
density_and_abdomen_centered = density_and_abdomen - np.mean(density_and_abdomen, axis = 0)
density_and_abdomen_centered.head(5)
Now when we do the approximation, things work much better.
density_and_abdomen_centered_rank_1_approximation = compute_rank_k_approximation(density_and_abdomen_centered, 1)
sns.scatterplot(x="density", y="abdomen", data=density_and_abdomen_centered)
sns.scatterplot(x="density", y="abdomen", data=density_and_abdomen_centered_rank_1_approximation);
Interpreting Principal Components¶
Let's revisit our child mortality and maternal fertility data from before.
sns.scatterplot(data = child_data, x = "mortality", y= "fertility")
plt.xlim([0, 14])
plt.ylim([0, 14])
plt.xticks(np.arange(0, 14, 2))
plt.yticks(np.arange(0, 14, 2));
Since we're going to be doing SVD, let's make sure to center our data first.
np.mean(child_data, axis = 0)
child_means = np.mean(child_data, axis = 0)
child_data_centered = child_data - child_means
sns.scatterplot(data = child_data_centered, x = "mortality", y= "fertility")
plt.xlim([-3, 11])
plt.ylim([-3, 11])
plt.xticks(np.arange(-3, 11, 2))
plt.yticks(np.arange(-3, 11, 2));
# plt.gcf().savefig("mortality_fertility_centered.png", dpi=300, bbox_inches="tight")
Principal Components and Singular Values¶
Singular Values and Variance¶
rectangle = pd.read_csv('data/rectangle_data.csv')
rectangle_centered = rectangle - np.mean(rectangle, axis = 0)
np.var(rectangle_centered)
sum(np.var(rectangle_centered))
u, s, vt = np.linalg.svd(rectangle_centered, full_matrices = False)
u[0:5, :]
np.diag(s)
s**2/rectangle_centered.shape[0]
sum(s**2/rectangle_centered.shape[0])
Practical PCA¶
Let's now step back and try to use PCA on our body measurement and congressional voting datasets.
body_data.head(5)
u, s, vt = np.linalg.svd(body_data, full_matrices = False)
We see that some of our principal components capture more variance than others.
s
Or we can compute the fraction of the variance captured by each principal component. The result seems shocking at first, as our data appears to be effectively rank 1.
np.round(s**2 / sum(s**2), 2)
This seems absurd, as clearly there are several variables that we expect to be show significant variation independent of each other, e.g. weight, height, and age. If this happens to you, it's probably because you forgot to center your data!
body_data_centered = body_data - np.mean(body_data, axis = 0)
body_data_centered.head(5)
u, s, vt = np.linalg.svd(body_data_centered, full_matrices = False)
This time, we see that the top singular value is no longer as dominant.
s
Looking now at the fraction of the variance captured by each principal component, we see that the top 2 or 3 principal components capture quite a lot of the variance.
np.round(s**2 / sum(s**2), 2)
We can also show this in the form of what is usually called a "scree plot".
plt.plot(s**2);
Thus, we expect that if we were to do a rank 3 approximation, we should get back data that's pretty close to where we started, as just those 3 dimensions capture 97% of the variance.
body_data_rank_3_approximation = compute_rank_k_approximation(body_data_centered, 3) + np.mean(body_data, axis = 0)
body_data_rank_3_approximation.head(5)
body_data.head(5)
One very interesting thing we can do is try to plot the principal components themselves. In this case, let's plot only the first two.
u, s, vt = np.linalg.svd(body_data_centered, full_matrices = False)
pcs = u * s
sns.scatterplot(x=pcs[:, 0], y=pcs[:, 1]);
np.argmax(pcs[:, 0])
body_data.iloc[38, :]
Election Data¶
from ds100_utils import fetch_and_cache
import yaml
from datetime import datetime
base_url = 'https://github.com/unitedstates/congress-legislators/raw/master/'
legislators_path = 'legislators-current.yaml'
f = fetch_and_cache(base_url + legislators_path, legislators_path)
legislators_data = yaml.load(open(f))
def to_date(s):
return datetime.strptime(s, '%Y-%m-%d')
legs = pd.DataFrame(
columns=['leg_id', 'first', 'last', 'gender', 'state', 'chamber', 'party', 'birthday'],
data=[[x['id']['bioguide'],
x['name']['first'],
x['name']['last'],
x['bio']['gender'],
x['terms'][-1]['state'],
x['terms'][-1]['type'],
x['terms'][-1]['party'],
to_date(x['bio']['birthday'])] for x in legislators_data])
legs.head(3)
# February 2019 House of Representatives roll call votes
# Downloaded using https://github.com/eyeseast/propublica-congress
votes = pd.read_csv('data/votes.csv')
votes = votes.astype({"roll call": str})
votes.head()
def was_yes(s):
if s.iloc[0] == 'Yes':
return 1
else:
return 0
vote_pivot = votes.pivot_table(index='member',
columns='roll call',
values='vote',
aggfunc=was_yes,
fill_value=0)
print(vote_pivot.shape)
vote_pivot.head()
vote_pivot_centered = vote_pivot - np.mean(vote_pivot, axis = 0)
vote_pivot_centered.head(5)
u, s, vt = np.linalg.svd(vote_pivot_centered, full_matrices = False)
np.round(s**2 / sum(s**2), 2)
plt.plot(s**2)
pcs = u * s
sns.scatterplot(x=pcs[:, 0], y=pcs[:, 1]);
vote2d = pd.DataFrame({
'member': vote_pivot.index,
'pc1': pcs[:, 0],
'pc2': pcs[:, 1]
}).merge(legs, left_on='member', right_on='leg_id')
vote2d[vote2d['pc1'] < 0]['party'].value_counts()
#top right only
vote2d.query('pc2 > -2 and pc1 > 0')['party'].value_counts()
sns.scatterplot(x="pc1", y="pc2", hue="party", data = vote2d);
vote2d['pc1_jittered'] = vote2d['pc1'] + np.random.normal(loc = 0, scale = 0.1, size = vote2d.shape[0])
vote2d['pc2_jittered'] = vote2d['pc2'] + np.random.normal(loc = 0, scale = 0.1, size = vote2d.shape[0])
sns.scatterplot(x="pc1_jittered", y="pc2_jittered", hue="party", data = vote2d);
vote2d.head(5)
vote2d[vote2d['pc1'] > 0]['party'].value_counts()
vote2d[vote2d['pc2'] < -1]
df = votes[votes['member'].isin(vote2d[vote2d['pc2'] < -1]['member'])]
df.groupby(['member', 'vote']).size()
legs.query("leg_id == 'A000367'")
Let's look at only people who have voted more than 15 times lately¶
votes.head(5)
num_yes_or_no_votes_per_member = votes.query("vote == 'Yes' or vote == 'No'").groupby("member").size()
num_yes_or_no_votes_per_member.head(5)
vote_pivot_with_yes_no_count = vote_pivot.merge(num_yes_or_no_votes_per_member.to_frame(), left_index = True, right_index = True, how="outer", ).fillna(0)
vote_pivot_with_yes_no_count = vote_pivot_with_yes_no_count.rename(columns = {0: 'yes_no_count'})
vote_pivot_with_yes_no_count.head(5)
regulars = vote_pivot_with_yes_no_count.query('yes_no_count >= 30')
regulars = regulars.drop('yes_no_count', 1)
regulars.shape
regulars_centered = regulars - np.mean(regulars, axis = 0)
regulars_centered.head(5)
u, s, vt = np.linalg.svd(regulars_centered, full_matrices = False)
np.round(s**2 / sum(s**2), 2)
pcs = u * s
sns.scatterplot(x=pcs[:, 0], y=pcs[:, 1]);
vote2d = pd.DataFrame({
'member': regulars_centered.index,
'pc1': pcs[:, 0],
'pc2': pcs[:, 1]
}).merge(legs, left_on='member', right_on='leg_id')
vote2d['pc1_jittered'] = vote2d['pc1'] + np.random.normal(loc = 0, scale = 0.1, size = vote2d.shape[0])
vote2d['pc2_jittered'] = vote2d['pc2'] + np.random.normal(loc = 0, scale = 0.1, size = vote2d.shape[0])
sns.scatterplot(x="pc1_jittered", y="pc2_jittered", hue="party", data = vote2d);
Exploring $V^T$¶
We can also look at Vt directly to try to gain insight into why each component is as it is.
num_votes = vt.shape[1]
votes = regulars.columns
def plot_pc(k):
plt.bar(votes, vt[k, :], alpha=0.7)
plt.xticks(votes, rotation=90);
with plt.rc_context({"figure.figsize": (12, 4)}):
plot_pc(0)
# plot_pc(1)
def compute_rank_k_approximation(data, k):
u, s, vt = np.linalg.svd(data, full_matrices = False)
return pd.DataFrame(u[:, 0:k] @ np.diag(s[0:k]) @ vt[0:k, :], columns = data.columns)
PCA vs. Linear Regression (Bonus)¶
Now we'll finally turn to trying to understand what the principle components and the low rank approximations actually mean!
Returning to our child mortality data from before, if we zero-center the child data, we see get back a better result. Note that we have to add back in the mean of each column to get things back into the right units.
means = np.mean(child_data, axis = 0)
child_data_centered = child_data - np.mean(child_data, axis = 0)
child_data_rank_1_approximation = compute_rank_k_approximation(child_data_centered, 1) + means
sns.scatterplot('mortality', 'fertility', data=child_data)
sns.scatterplot('mortality', 'fertility', data=child_data_rank_1_approximation)
plt.legend(['data', 'rank 1 approx'])
plt.xlim([0, 14])
plt.ylim([0, 14])
plt.xticks(np.arange(0, 14, 2))
plt.yticks(np.arange(0, 14, 2));
We can also give our rank 1 approximation as a line, showing the 1D subspace (in black) that our data is being projected onto.
sns.scatterplot('mortality', 'fertility', data=child_data)
sns.lineplot('mortality', 'fertility', data=child_data_rank_1_approximation, color='black')
plt.legend(['rank 1 approx', 'data']);
plt.xlim([0, 14])
plt.ylim([0, 14])
plt.xticks(np.arange(0, 14, 2))
plt.yticks(np.arange(0, 14, 2));
This plot probably brings to mind linear regression from Data 8. But PCA is NOT the same thing linear regression. Let's plot the regression lines for this data for comparison. Recall that the regression line gives us, e.g. the best possible linear prediction of the fertility given the mortality.
x, y = child_data['mortality'], child_data['fertility']
slope_x, intercept_x = np.polyfit(x, y, 1) # simple linear regression
scatter14(child_data)
plt.plot(x, slope_x * x + intercept_x, c = 'g')
for _, row in child_data.sample(20).iterrows():
tx, ty = row['mortality'], row['fertility']
plt.plot([tx, tx], [slope_x * tx + intercept_x, ty], c='red')
plt.legend(['predicted fertility']);
In the plot above, the green regression line given minimizes the sum of the squared errors, given as red vertical lines.
We could also do the opposite thing, and try to predict fertility from mortality.
x, y = child_data['mortality'], child_data['fertility']
slope_y, intercept_y = np.polyfit(y, x, 1) # simple linear regression
scatter14(child_data)
plt.plot(slope_y * y + intercept_y, y, c = 'purple')
for _, row in child_data.sample(20).iterrows():
tx, ty = row['mortality'], row['fertility']
plt.plot([tx, slope_y * ty + intercept_y], [ty, ty], c='red')
plt.legend(['predicted mortality']);
In the plot above, the green regression line given minimizes the sum of the squared errors, given as red horizontal lines.
Plotting the two regression lines and the 1D subspace chosen by PCA, we see that all 3 are distinct!
sns.lineplot('mortality', 'fertility', data=child_data_rank_1_approximation, color="black")
plt.plot(x, slope_x * x + intercept_x, c = 'g')
plt.plot(slope_y * y + intercept_y, y, c = 'purple');
sns.scatterplot('mortality', 'fertility', data=child_data)
plt.legend(['rank 1 approx', 'predicted fertility', 'predicted mortality'])
plt.xlim([0, 14])
plt.ylim([0, 14])
plt.xticks(np.arange(0, 14, 2))
plt.yticks(np.arange(0, 14, 2));
Given that the green line minimizes the "vertical error" and the purple line minimizes the "horizontal error". You might wonder what the black line minimizes. It turns out, it minimizes the "diagonal" error, i.e. the error in the direction perpendicular to itself.
sns.lineplot('mortality', 'fertility', data=child_data_rank_1_approximation, color="black")
plt.plot(x, slope_x * x + intercept_x, c = 'g')
plt.plot(slope_y * y + intercept_y, y, c = 'purple');
sns.scatterplot('mortality', 'fertility', data=child_data)
for idx, tdata in child_data.reset_index().sample(20).iterrows():
tx = tdata["mortality"]
ty = tdata["fertility"]
tx_projected = child_data_rank_1_approximation.iloc[idx, 0]
ty_projected = child_data_rank_1_approximation.iloc[idx, 1]
plt.plot([tx, tx_projected], [ty, ty_projected], c='red')
plt.xlim([0, 14])
plt.ylim([0, 14])
plt.xticks(np.arange(0, 14, 2))
plt.yticks(np.arange(0, 14, 2));
plt.legend(['rank 1 approx', 'predicted fertility', 'predicted mortality']);
The function in the following cell makes it easy to make similar plots for whatever dataset you might be interested in.
def plot_x_regression_y_regression_1d_approximation(data):
xname = data.columns[0]
yname = data.columns[1]
x, y = data[xname], data[yname]
slope_x, intercept_x = np.polyfit(x, y, 1) # simple linear regression
x, y = data[xname], data[yname]
slope_y, intercept_y = np.polyfit(y, x, 1) # simple linear regression
means = np.mean(data, axis = 0)
rank_1_approximation = compute_rank_k_approximation(data - means, 1) + means
sns.lineplot(x=xname, y=yname, data=rank_1_approximation, color="black")
plt.plot(x, slope_x * x + intercept_x, c = 'g')
plt.plot(slope_y * y + intercept_y, y, c = 'purple');
sns.scatterplot(xname, yname, data=data)
for idx, tdata in data.reset_index().sample(20).iterrows():
tx = tdata[xname]
ty = tdata[yname]
tx_projected = rank_1_approximation.iloc[idx, 0]
ty_projected = rank_1_approximation.iloc[idx, 1]
plt.plot([tx, tx_projected], [ty, ty_projected], c='red')
plt.legend(['1D PCA Subspace', 'predicted ' + xname, 'predicted ' + yname])
plot_x_regression_y_regression_1d_approximation(body_data.drop(41)[["adiposity", "bicep"]])
