Live Lecture 4 Supplemental Notebook¶
Data 100, Summer 2020¶
Suraj Rampure and Allen Shen
This notebook has 5 sections:
- An overview of the modeling process, and how it parallels with
sklearn - An exploration of how training RMSE doesn't decrease when we add more features, but test RMSE can
- One-hot encoding, and issues with an intercept term
- Redundant features and rank
- What makes models linear
In [1]:
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
import sklearn.linear_model as lm
1. An overview of the modeling process¶
- Choose a model
- Choose a loss function
model = lm.LinearRegression()does both of these for us!
- Minimize average loss to find the optimal $\hat{\theta}$
model.fit(X, y)model.coef_andmodel.intercept_give us the values of $\hat{\theta}$ after fitting
Use our model to make predictions: $\hat{\mathbb{Y}} = \mathbb{X} \hat{\theta}$
model.predict(X)
2. RMSE on training data never increases by adding more features¶
In [2]:
df = sns.load_dataset('tips')
df.head()
Out[2]:
In [3]:
def rmse(y, yhat):
return np.sqrt(np.mean((y - yhat)**2))
Let's start by fitting a simple linear regression model to our familiar tips data. Specifically, we will use total_bill to predict tip.
In [4]:
model1 = lm.LinearRegression()
model1.fit(df[['total_bill']], df['tip'])
pred1 = model1.predict(df[['total_bill']])
In [5]:
rmse(df['tip'], pred1)
Out[5]:
In [6]:
model1.coef_
Out[6]:
In [7]:
model1.intercept_
Out[7]:
In [8]:
plt.scatter(df['tip'], pred1);
Notably, our RMSE was 1.0178504025697377.
Now, let's add a completely unrelated column to our data, and include it as a feature in our model.
In [9]:
df['useless'] = np.random.randn(len(df)) * 342
In [10]:
df
Out[10]:
In [11]:
model2 = lm.LinearRegression()
model2.fit(df[['total_bill', 'useless']], df['tip'])
pred2 = model2.predict(df[['total_bill', 'useless']])
In [12]:
rmse(df['tip'], pred2)
Out[12]:
In [13]:
model2.coef_
Out[13]:
In [14]:
model2.intercept_
Out[14]:
In [15]:
plt.scatter(df['tip'], pred2);
Our new RMSE was marginally lower! Why?
Note that the coefficient for our useless feature is very close to 0.
What about Multiple $R^2$?¶
For the original model:
In [16]:
np.corrcoef(pred1, df['tip'])[0, 1]**2
Out[16]:
In [17]:
np.var(pred1) / np.var(df['tip'])
Out[17]:
Note: model.score for a LinearRegression model also computes the $R^2$ value! See:
In [18]:
model1.score(df[['total_bill']], df['tip'])
Out[18]:
You should note that the above three are all the same.
For the model with the useless feature:
In [19]:
np.var(pred2) / np.var(df['tip'])
Out[19]:
Recall, we can interpret $R^2$ as being the proportion of variance in our true $y$ values that our fitted values capture. This is saying that our model with the useless feature included accounts for more of the variation in our true $y$ values than the first model does.
Does this make sense?
Let's look at how such a model performs on unseen ("test") data.¶
We haven't yet formally taught you how to use scikit-learn's inbuilt train/test split, so we will do this by hand for now.
In [20]:
idx = np.arange(len(df))
np.random.shuffle(idx)
In [21]:
idx
Out[21]:
In [22]:
len(idx)
Out[22]:
In [23]:
split_point = int((3/4) * len(idx))
In [24]:
split_point
Out[24]:
In [25]:
train, test = df.iloc[idx[:split_point]], df.iloc[idx[split_point:]]
In [26]:
train
Out[26]:
In [27]:
test
Out[27]:
In [28]:
new_model1 = lm.LinearRegression()
new_model1.fit(train[['total_bill']], train['tip'])
new_pred1_train = new_model1.predict(train[['total_bill']])
new_pred1_test = new_model1.predict(test[['total_bill']])
In [29]:
new_model1_train_rmse = rmse(train['tip'], new_pred1_train)
new_model1_test_rmse = rmse(test['tip'], new_pred1_test)
new_model1_train_rmse, new_model1_test_rmse
Out[29]:
Now, for our model with the useless feature:
In [30]:
new_model2 = lm.LinearRegression()
new_model2.fit(train[['total_bill', 'useless']], train['tip'])
new_pred2_train = new_model2.predict(train[['total_bill', 'useless']])
new_pred2_test = new_model2.predict(test[['total_bill', 'useless']])
In [31]:
new_model2_train_rmse = rmse(train['tip'], new_pred2_train)
new_model2_test_rmse = rmse(test['tip'], new_pred2_test)
new_model2_train_rmse, new_model2_test_rmse
Out[31]:
In [32]:
new_model1.coef_
Out[32]:
In [33]:
new_model2.coef_
Out[33]:
Note, here training RMSE went down, but test RMSE went up. This is generally what happens when you include features that aren't truly relevant to the underlying relationship in your data. We call this overfitting.
3. One hot encoding¶
In [34]:
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
import sklearn.linear_model as lm
from sklearn.feature_extraction import DictVectorizer
In [35]:
df = sns.load_dataset('tips')
df.head()
Out[35]:
Why do we need one hot encoding?¶
In [36]:
model1 = lm.LinearRegression()
model1.fit(df.drop(columns='tip'), df['tip'])
How to perform a one hot encoding in scikit-learn¶
In [37]:
df = df.drop(columns='tip')
cat_cols = ['sex', 'smoker', 'day', 'time']
In [38]:
vec_enc = DictVectorizer()
vec_enc.fit(df[cat_cols].to_dict(orient='records'))
cat_data = vec_enc.transform(df[cat_cols].to_dict(orient='records')).toarray()
cat_data
Out[38]:
In [39]:
cat_data_names = vec_enc.get_feature_names()
cat_data_names
Out[39]:
In [40]:
cat_data = pd.DataFrame(cat_data, columns=cat_data_names)
df_ohe = pd.concat([df, cat_data], axis=1).drop(columns=cat_cols) # Drop original categorical columns
df_ohe.head()
Out[40]:
In [41]:
df_ohe.shape
Out[41]:
In [42]:
X = pd.concat([df_ohe, pd.Series(np.ones(df_ohe.shape[0]), name='intercept')], axis=1)
X.head()
Out[42]:
In [43]:
X.shape
Out[43]:
What does this output?¶
In [44]:
np.linalg.matrix_rank(X)
Out[44]:
What's the issue with the design matrix above?¶
In [45]:
X = X.drop(columns=['day=Sat', 'sex=Female', 'smoker=No', 'time=Dinner'])
X.head()
Out[45]:
In [46]:
X.shape
Out[46]:
What does this output?¶
In [47]:
np.linalg.matrix_rank(X)
Out[47]:
One hot encoding with Pandas¶
In [48]:
df_ohe2 = pd.get_dummies(df)
df_ohe2.head()
Out[48]:
In [49]:
X_2 = pd.concat([df_ohe2, pd.Series(np.ones(df_ohe2.shape[0]), name='intercept')], axis=1)
X_2.head()
Out[49]:
In [50]:
X_2.shape
Out[50]:
In [51]:
np.linalg.matrix_rank(X_2)
Out[51]:
In [52]:
df_ohe2 = pd.get_dummies(df, drop_first=True)
df_ohe2.head()
Out[52]:
In [53]:
X_2 = pd.concat([df_ohe2, pd.Series(np.ones(df_ohe2.shape[0]), name='intercept')], axis=1)
X_2.head()
Out[53]:
In [54]:
X_2.shape
Out[54]:
In [55]:
np.linalg.matrix_rank(X_2)
Out[55]:
More examples of one hot encoding¶
In [56]:
X_3 = pd.get_dummies(df[['total_bill', 'sex']])
X_3.head()
Out[56]:
What would this output?¶
In [57]:
np.linalg.matrix_rank(X_3)
Out[57]:
In [58]:
X_4 = pd.concat([X_3, pd.Series(np.ones(X_3.shape[0]), name='intercept')], axis=1)
X_4.head()
Out[58]:
What would this output?¶
In [59]:
np.linalg.matrix_rank(X_4)
Out[59]:
In [60]:
(X_4['intercept'] - X_4['sex_Male']).iloc[:5]
Out[60]:
In [61]:
X_5 = pd.get_dummies(df[['total_bill', 'sex', 'smoker']])
X_5.tail(5)
Out[61]:
What would this output?¶
In [62]:
np.linalg.matrix_rank(X_5)
Out[62]:
In [63]:
(X_5['sex_Male'] + X_5['sex_Female'] - X_5['smoker_Yes']).iloc[-5:]
Out[63]:
In [64]:
X_6 = pd.concat([X_5, pd.Series(np.ones(X_5.shape[0]), name='intercept')], axis=1)
X_6.tail()
Out[64]:
What would this output?¶
In [65]:
np.linalg.matrix_rank(X_6)
Out[65]:
4. Duplicate features¶
In [66]:
df = sns.load_dataset('tips')
In [67]:
X_7 = df[['total_bill', 'size', 'size']]
X_7 = pd.concat([X_7, pd.Series(np.ones(X_7.shape[0]), name='intercept')], axis=1)
X_7.head()
Out[67]:
What would this output?¶
In [68]:
np.linalg.matrix_rank(X_7)
Out[68]:
What's the issue with this again?¶
In [69]:
model2 = lm.LinearRegression(fit_intercept=False)
model2.fit(X_7, df['tip'])
Out[69]:
In [70]:
model2_coef = model2.coef_
model2_coef
Out[70]:
In [71]:
model2.intercept_
Out[71]:
In [72]:
model2.predict(X_7)[:5]
Out[72]:
In [73]:
X_7.iloc[:5] @ model2_coef
Out[73]:
How can I change the model coefficients so that the predictions are the same?¶
Our model is:
$$\theta_0 + \theta_1 \cdot size + \theta_2 \cdot size + \theta_3 \cdot total\_bill$$In [74]:
model2_coef_modified = model2_coef.copy()
model2_coef_modified[1] = model2_coef[1] - 1000000000
model2_coef_modified[2] = model2_coef[2] + 1000000000
In [75]:
model2_coef
Out[75]:
In [76]:
model2_coef_modified
Out[76]:
In [77]:
X_7.iloc[:5] @ model2_coef_modified
Out[77]:
Our model is now:
$$\theta_0 + (\theta_1 - 100) \cdot size + (\theta_2 + 100) \cdot size + \theta_3 \cdot total\_bill$$$$= \theta_0 + \theta_1 \cdot size - 100 \cdot size + \theta_2 * size + 100 \cdot size + \theta_3 \cdot total\_bill$$$$= \theta_0 + \theta_1 \cdot size + \theta_2 \cdot size + \theta_3 \cdot total\_bill$$which is the same as before!
Using 2 times size as a feature¶
In [78]:
X_8 = df[['total_bill', 'size']]
X_8.loc[:, '2 * size'] = 2 * X_8['size']
X_8 = pd.concat([X_8, pd.Series(np.ones(X_8.shape[0]), name='intercept')], axis=1)
X_8.head()
Out[78]:
What would this output?¶
In [79]:
np.linalg.matrix_rank(X_8)
Out[79]:
In [80]:
model3 = lm.LinearRegression(fit_intercept=False)
model3.fit(X_8, df['tip'])
Out[80]:
In [81]:
model3_coef = model3.coef_
model3_coef
Out[81]:
In [82]:
model3.predict(X_8)[:5]
Out[82]:
In [83]:
X_8.iloc[:5] @ model3_coef
Out[83]:
How can I change the coefficients so that the predictions are the same?¶
Can I do the same thing as before?
In [84]:
model3_coef_modified = model3_coef.copy()
model3_coef_modified[1] = model3_coef[1] - 100
model3_coef_modified[2] = model3_coef[2] + 100
In [85]:
X_8.iloc[:5] @ model3_coef_modified
Out[85]:
In [86]:
model3_coef_modified = model3_coef.copy()
model3_coef_modified[1] = model3_coef[1] - 100
model3_coef_modified[2] = model3_coef[2] + 50
In [87]:
X_8.iloc[:5] @ model3_coef_modified
Out[87]:
Thought exercise: What happens if I try to add 2 * size + 3 as a feature?¶
In [88]:
X_9 = df[['total_bill', 'size']]
X_9.loc[:, '2 * size + 3'] = 2 * X_9['size'] + 3
X_9 = pd.concat([X_9, pd.Series(np.ones(X_8.shape[0]), name='intercept')], axis=1)
X_9.head()
Out[88]:
What would this output?¶
In [89]:
np.linalg.matrix_rank(X_9)
Out[89]:
Adding size squared as a feature¶
In [90]:
X_10 = df[['total_bill', 'size']]
X_10.loc[:, 'size ** 2'] = X_10['size'] ** 2
X_10 = pd.concat([X_10, pd.Series(np.ones(X_10.shape[0]), name='intercept')], axis=1)
X_10.head()
Out[90]:
What would this output?¶
In [91]:
np.linalg.matrix_rank(X_10)
Out[91]:
5. What makes a model linear?¶
Is the following model linear? (Suppose $x$ represents a single observation of our raw data matrix.)
$$f_\theta(x) = \theta_0 + \theta_1 x_1 + \theta_2 x_2 + \theta_3 \sin(x_1) + \theta_4 \cos(x_1x_2) e^{x_1}$$Yes, because it is linear in terms of the parameters. We could formulate this model as $$f_\theta(x) = x^T \theta$$ where $x = \begin{bmatrix} 1 \\ x_1 \\ x_2 \\ \sin(x_1) \\ \cos(x_1x_2) e^{x_1} \end{bmatrix}$ and $\theta = \begin{bmatrix} \theta_0 \\ \theta_1 \\ \theta_2 \\ \theta_3 \\ \theta_4 \end{bmatrix}$.
What about this model?
$$f_\theta(x) = \theta_0 + \theta_1 x + \theta_2 \sin(\theta_3 x)$$No, because $\theta_3$ is within a $\sin$ function. We cannot write this model in the form $x^T \theta$, so it is not a linear model. It is still a model, though, and we can find its optimal parameters, but just not using least squares.