In [1]:

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn import tree, model_selection, datasets, metrics, ensemble, linear_model, neighbors
import graphviz as gv
import seaborn as sns

In the lab, we applied random forests to the Boston data using mtry=6 and using ntree=25 and ntree=500. Create a plot displaying the test error resulting from random forests on this data set for a more comprehensive range of values for mtry and ntree. You can model your plot after Figure 8.10. Describe the results obtained.

In [2]:

boston_df = datasets.load_boston()
boston_df = pd.DataFrame(data=np.c_[boston_df['data'], boston_df['target']], columns= [c for c in boston_df['feature_names']] + ['Price'])
boston_df.head()

Out[2]:

	CRIM	ZN	INDUS	NOX	RM	AGE	DIS	RAD	TAX	PTRATIO	B	LSTAT	Price
0	0.00632	18.0	2.31	0.538	6.575	65.2	4.0900	1.0	296.0	15.3	396.90	4.98	24.0
1	0.02731	0.0	7.07	0.469	6.421	78.9	4.9671	2.0	242.0	17.8	396.90	9.14	21.6
2	0.02729	0.0	7.07	0.469	7.185	61.1	4.9671	2.0	242.0	17.8	392.83	4.03	34.7
3	0.03237	0.0	2.18	0.458	6.998	45.8	6.0622	3.0	222.0	18.7	394.63	2.94	33.4
4	0.06905	0.0	2.18	0.458	7.147	54.2	6.0622	3.0	222.0	18.7	396.90	5.33	36.2

In [2]:

def run_cv(df, n_splits, extract_x, extract_y, fit_model, loss_fn):
    cv_result = pd.DataFrame()
    for train_idx, test_idx in model_selection.KFold(n_splits=n_splits).split(df):
        train, test = df.iloc[train_idx], df.iloc[test_idx]
        errors = []
        # train
        train_X = extract_x(train)
        train_y = extract_y(train)
        model = fit_model(train_X, train_y)
        # test
        test_X =  extract_x(test)
        test_y = extract_y(test)
        preds = model.predict(test_X)
        errors.append(loss_fn(preds, test_y))
        cv_result = cv_result.append(pd.Series(errors), ignore_index=True)
    return cv_result

In [78]:

rnge = range(1, 50)

sqrt_rmses = [run_cv(boston_df,
                     5,
                     lambda df: df.drop('Price', axis=1),
                     lambda df:  df.Price,
                     lambda x, y: ensemble.RandomForestRegressor(n_estimators=i, max_features='sqrt').fit(x,y),
                     lambda preds, true: np.sqrt(metrics.mean_squared_error(true, preds))).mean()[0] for i in rnge]

In [ ]:

log2_rmses = [run_cv(boston_df,
                     5,
                     lambda df: df.drop('Price', axis=1),
                     lambda df:  df.Price,
                     lambda x, y: ensemble.RandomForestRegressor(n_estimators=i, max_features='log2').fit(x,y),
                     lambda preds, true: np.sqrt(metrics.mean_squared_error(true, preds))).mean()[0] for i in rnge]

In [76]:

all_fx_rmses = [run_cv(boston_df,
                     5,
                     lambda df: df.drop('Price', axis=1),
                     lambda df:  df.Price,
                     lambda x, y: ensemble.RandomForestRegressor(n_estimators=i, max_features=None).fit(x,y),
                     lambda preds, true: np.sqrt(metrics.mean_squared_error(true, preds))).mean()[0] for i in rnge]

In [77]:

df = pd.DataFrame({'sqrt_rmses' : sqrt_rmses,
                   'log2_rmses' : log2_rmses,
                   'all_fx_rmses': all_fx_rmses,
                   'trees' : rnge})
sns.lineplot(x='trees', y='sqrt_rmses', data=df, color='tab:blue')
sns.lineplot(x='trees', y='log2_rmses', data=df, color='tab:green')
sns.lineplot(x='trees', y='all_fx_rmses', data=df, color='tab:red')

Out[77]:

<matplotlib.axes._subplots.AxesSubplot at 0x1128a1278>

In the lab, a classification tree was applied to the Carseats data set after converting Sales into a qualitative response variable. Now we will seek to predict Sales using regression trees and related approaches, treating the response as a quantitative variable.

In [116]:

car_df = pd.read_csv('carseats.csv')
car_df = car_df.drop(car_df.columns[0], axis=1)
car_df['Urban'] = car_df['Urban'] == 'Yes'
car_df['US'] = car_df['US'] == 'Yes'
car_df['ShelveLoc'] = car_df['ShelveLoc'].map({'Bad' : 0, 'Medium': 1, 'Good' : 2})
car_df.head()

Out[116]:

	Sales	CompPrice	Income	Advertising	Population	Price	ShelveLoc	Age	Education	Urban	US
0	9.50	138	73	11	276	120	0	42	17	True	True
1	11.22	111	48	16	260	83	2	65	10	True	True
2	10.06	113	35	10	269	80	1	59	12	True	True
3	7.40	117	100	4	466	97	1	55	14	True	True
4	4.15	141	64	3	340	128	0	38	13	True	False

(a) Split the data set into a training set and a test set.

In [118]:

train, test = model_selection.train_test_split(car_df, test_size=0.2)

(b) Fit a regression tree to the training set. Plot the tree, and interpret the results. What test MSE do you obtain?

In [119]:

reg = tree.DecisionTreeRegressor(max_depth=3) 
train_x = train.drop('Sales', axis=1)
train_y = train.Sales
reg.fit(train_x,train_y)
test_x = test.drop('Sales', axis=1)
test_y = test.Sales
preds = reg.predict(test_x)
rmse = np.sqrt(metrics.mean_squared_error(test_y, preds))
print('test rmse: {}'.format(rmse))

test rmse: 2.2041731241021694

In [80]:

dot_dat = tree.export_graphviz(reg, out_file=None, rotate=True)
graph = gv.Source(dot_dat)
graph

Out[80]:

(c) Use cross-validation in order to determine the optimal level of tree complexity. Does pruning the tree improve the test MSE?

.. omitting this as tree pruning isn't easily available in python world, the python community prefer to control variance with boosting, bagging and random forest.

(d) Use the bagging approach in order to analyze this data. What test MSE do you obtain? Use the importance() function to de- termine which variables are most important.

In [83]:

reg = ensemble.RandomForestRegressor(n_estimators=100, max_features=None)
train_x = train.drop('Sales', axis=1)
train_y = train.Sales
reg.fit(train_x,train_y)
test_x = test.drop('Sales', axis=1)
test_y = test.Sales
preds = reg.predict(test_x)
rmse = np.sqrt(metrics.mean_squared_error(test_y, preds))
print('training rmse: {}'.format(rmse))

_,_ = plt.subplots(figsize=(15, 7.5))
bar_df = pd.DataFrame({'predictor': train_x.columns, 'importance' : reg.feature_importances_})
sns.barplot(x='predictor', y='importance', data=bar_df.sort_values('importance', ascending=False))

training rmse: 1.5726550643736208

Out[83]:

<matplotlib.axes._subplots.AxesSubplot at 0x11289b940>

(e) Use random forests to analyze this data. What test MSE do you obtain? Use the importance() function to determine which variables are most important. Describe the effect of m, the number of variables considered at each split, on the error rate obtained.

In [85]:

reg = ensemble.RandomForestRegressor(n_estimators=100, max_features='sqrt')
train_x = train.drop('Sales', axis=1)
train_y = train.Sales
reg.fit(train_x,train_y)
test_x = test.drop('Sales', axis=1)
test_y = test.Sales
preds = reg.predict(test_x)
rmse = np.sqrt(metrics.mean_squared_error(test_y, preds))
print('training rmse: {}'.format(rmse))

_,_ = plt.subplots(figsize=(15, 7.5))
bar_df = pd.DataFrame({'predictor': train_x.columns, 'importance' : reg.feature_importances_})
sns.barplot(x='predictor', y='importance', data=bar_df.sort_values('importance', ascending=False))

training rmse: 1.7675543885549883

Out[85]:

<matplotlib.axes._subplots.AxesSubplot at 0x112b8ca58>

In [103]:

rnge = np.linspace(0.1, 1, 10)

sqrt_rmses = [run_cv(car_df,
                     10,
                     lambda df: df.drop('Sales', axis=1),
                     lambda df:  df.Sales,
                     lambda x, y: ensemble.RandomForestRegressor(n_estimators=100, max_features=i).fit(x,y),
                     lambda preds, true: np.sqrt(metrics.mean_squared_error(true, preds))).mean()[0] for i in rnge]

In [104]:

line_df = pd.DataFrame({'m' : rnge * car_df.shape[1] - 1, 'rmse' : sqrt_rmses})
sns.lineplot(x='m', y='rmse', data=line_df)

Out[104]:

<matplotlib.axes._subplots.AxesSubplot at 0x1132a1940>

This problem involves the OJ data set which is part of the ISLR package.

In [6]:

oj_df = pd.read_csv('oj.csv')
oj_df = oj_df.drop(oj_df.columns[0], axis=1)
oj_df.Purchase = oj_df.Purchase.map({'CH' : 1, 'MM': 0})
oj_df.Store7 = oj_df.Store7.map({'Yes' : 1, 'No': 0})
oj_df.head()

Out[6]:

	Purchase	WeekofPurchase	StoreID	PriceCH	PriceMM	DiscCH	DiscMM	SpecialMM	LoyalCH	SalePriceMM	SalePriceCH	PriceDiff	Store7	PctDiscMM	PctDiscCH	ListPriceDiff	STORE
0	1	237	1	1.75	1.99	0.00	0.0	0	0.500000	1.99	1.75	0.24	0	0.000000	0.000000	0.24	1
1	1	239	1	1.75	1.99	0.00	0.3	1	0.600000	1.69	1.75	-0.06	0	0.150754	0.000000	0.24	1
2	1	245	1	1.86	2.09	0.17	0.0	0	0.680000	2.09	1.69	0.40	0	0.000000	0.091398	0.23	1
3	0	227	1	1.69	1.69	0.00	0.0	0	0.400000	1.69	1.69	0.00	0	0.000000	0.000000	0.00	1
4	1	228	7	1.69	1.69	0.00	0.0	0	0.956535	1.69	1.69	0.00	1	0.000000	0.000000	0.00	0

(a) Create a training set containing a random sample of 800 observations, and a test set containing the remaining observations.

In [133]:

train, test = model_selection.train_test_split(oj_df, test_size=oj_df.shape[0] - 800)

(b) Fit a tree to the training data, with Purchase as the response and the other variables as predictors. Use the summary() function to produce summary statistics about the tree, and describe the results obtained. What is the training error rate? How many terminal nodes does the tree have?

In [162]:

reg = tree.DecisionTreeClassifier(max_depth=3) 
train_x = train.drop('Purchase', axis=1)
train_y = train.Purchase
reg.fit(train_x,train_y)

# train error
preds = reg.predict(train_x)
train_error = 1 - ((train_y == preds).sum() / train_y.shape[0])
print('train error: {}'.format(train_error))

train error: 0.16249999999999998

(c) Type in the name of the tree object in order to get a detailed text output. Pick one of the terminal nodes, and interpret the information displayed. (d) Create a plot of the tree, and interpret the results.

In [138]:

dot_dat = tree.export_graphviz(reg, out_file=None, rotate=True)
graph = gv.Source(dot_dat)
graph

Out[138]:

(e) Predict the response on the test data, and produce a confusion matrix comparing the test labels to the predicted test labels. What is the test error rate?

In [163]:

test_x = test.drop('Purchase', axis=1)
test_y = test.Purchase
preds = reg.predict(test_x)
conf = pd.DataFrame(metrics.confusion_matrix(test_y, preds))
conf

Out[163]:

	0	1
0	75	22
1	27	146

In [164]:

print('test error rate: {}'.format( 1 - ((conf[0][0] + conf[1][1]) / test_y.shape[0])))

test error rate: 0.18148148148148147

(f) Apply the cv.tree() function to the training set in order to determine the optimal tree size.

.. not easily available in python

(g) Produce a plot with tree size on the x-axis and cross-validated classification error rate on the y-axis. (h) Which tree size corresponds to the lowest cross-validated classification error rate?

In [8]:

rnge = range(1, 50)

errors = [run_cv(oj_df,
                 5,
                 lambda df: df.drop('Purchase', axis=1),
                 lambda df:  df.Purchase,
                 lambda x, y: tree.DecisionTreeClassifier(max_depth=i) .fit(x,y),
                 lambda preds, true: 1 - ((preds == true).sum() / preds.shape[0])).mean()[0] for i in rnge]

In [14]:

line_df = pd.DataFrame({'max_depth' : rnge, 'error' : errors})
sns.lineplot(x='max_depth', y='error', data=line_df)

Out[14]:

<matplotlib.axes._subplots.AxesSubplot at 0x1162366d8>

(i) Produce a pruned tree corresponding to the optimal tree size obtained using cross-validation. If cross-validation does not lead to selection of a pruned tree, then create a pruned tree with five terminal nodes. (j) Compare the training error rates between the pruned and un- pruned trees. Which is higher? (k) Compare the test error rates between the pruned and unpruned trees. Which is higher?

.. pruning not easily available in python world

We now use boosting to predict Salary in the Hitters data set.

In [40]:

hit_df = pd.read_csv('hitters.csv')
hit_df = hit_df.drop(hit_df.columns[0], axis = 1)
hit_df = hit_df.dropna()
hit_df.Salary = np.log(hit_df.Salary)
hit_df.League = hit_df.League.map({'N':1, 'A':0})
hit_df.Division = hit_df.Division.map({'W':1, 'E':0})
hit_df.NewLeague = hit_df.NewLeague.map({'N':1, 'A':0})
hit_df.head()

Out[40]:

	AtBat	Hits	HmRun	Runs	RBI	Walks	Years	CAtBat	CHits	CHmRun	CRuns	CRBI	CWalks	League	Division	PutOuts	Assists	Errors	Salary	NewLeague
1	315	81	7	24	38	39	14	3449	835	69	321	414	375	1	1	632	43	10	6.163315	1
2	479	130	18	66	72	76	3	1624	457	63	224	266	263	0	1	880	82	14	6.173786	0
3	496	141	20	65	78	37	11	5628	1575	225	828	838	354	1	0	200	11	3	6.214608	1
4	321	87	10	39	42	30	2	396	101	12	48	46	33	1	0	805	40	4	4.516339	1
5	594	169	4	74	51	35	11	4408	1133	19	501	336	194	0	1	282	421	25	6.620073	0

(b) Create a training set consisting of the first 200 observations, and a test set consisting of the remaining observations.

In [43]:

train, test = model_selection.train_test_split(hit_df, test_size=hit_df.shape[0] - 200)

(c) Perform boosting on the training set with 1,000 trees for a range of values of the shrinkage parameter λ. Produce a plot with different shrinkage values on the x-axis and the corresponding training set MSE on the y-axis.

In [51]:

def hitters_boosting_test_error(lr):
    train_x = train.drop('Salary', axis=1)
    train_y = train.Salary
    reg = ensemble.GradientBoostingRegressor(n_estimators=1000, learning_rate=lr).fit(train_x, train_y)
    test_x = test.drop('Salary', axis=1)
    test_y = test.Salary
    preds = reg.predict(test_x)
    return np.sqrt(metrics.mean_squared_error(test_y, preds))

lambdas = np.linspace(0.01, 1, 100)
test_errors = [hitters_boosting_test_error(lr) for lr in lambdas]

In [52]:

def hitters_boosting_training_error(lr):
    train_x = train.drop('Salary', axis=1)
    train_y = train.Salary
    reg = ensemble.GradientBoostingRegressor(n_estimators=1000, learning_rate=lr).fit(train_x, train_y)
    preds = reg.predict(train_x)
    return np.sqrt(metrics.mean_squared_error(train_y, preds))

lambdas = np.linspace(0.01, 1, 100)
train_errors = [hitters_boosting_training_error(lr) for lr in lambdas]

In [58]:

line_df = pd.DataFrame({'test_error': test_errors,
                        'train_error': train_errors,
                        'lambda' : lambdas})
sns.lineplot(x='lambda', y='test_error', data=line_df, color='tab:blue')
sns.lineplot(x='lambda', y='train_error', data=line_df, color='tab:red')

Out[58]:

<matplotlib.axes._subplots.AxesSubplot at 0x1167bed68>

(e) Compare the test MSE of boosting to the test MSE that results from applying two of the regression approaches seen in Chapters 3 and 6.

In [80]:

train_x = train.drop('Salary', axis=1)
train_y = train.Salary
reg = linear_model.Ridge(alpha=10000).fit(train_x, train_y)
test_x = test.drop('Salary', axis=1)
test_y = test.Salary
preds = reg.predict(test_x)
np.sqrt(metrics.mean_squared_error(test_y, preds))

Out[80]:

0.5973063683221087

In [92]:

train_x = train.drop('Salary', axis=1)
train_y = train.Salary
reg = linear_model.Lasso(alpha=0.5).fit(train_x, train_y)
test_x = test.drop('Salary', axis=1)
test_y = test.Salary
preds = reg.predict(test_x)
np.sqrt(metrics.mean_squared_error(test_y, preds))

Out[92]:

0.5953299537143887

(f) Which variables appear to be the most important predictors in the boosted model?

In [97]:

train_x = train.drop('Salary', axis=1)
train_y = train.Salary
reg = ensemble.GradientBoostingRegressor(n_estimators=1000, learning_rate=0.2).fit(train_x, train_y)

# plot
_,_ = plt.subplots(figsize=(15, 7.5))
bar_df = pd.DataFrame({'predictor': train_x.columns, 'importance' : reg.feature_importances_})
sns.barplot(x='predictor', y='importance', data=bar_df.sort_values('importance', ascending=False))

Out[97]:

<matplotlib.axes._subplots.AxesSubplot at 0x1173d5518>

(g) Now apply bagging to the training set. What is the test set MSE for this approach?

In [100]:

train_x = train.drop('Salary', axis=1)
train_y = train.Salary
reg = ensemble.RandomForestRegressor(n_estimators=100, max_features=None).fit(train_x,train_y)
test_x = test.drop('Salary', axis=1)
test_y = test.Salary
preds = reg.predict(test_x)
np.sqrt(metrics.mean_squared_error(test_y, preds))

Out[100]:

0.4969384424303783

This question uses the Caravan data set.

In [107]:

van_df = pd.read_csv('caravan.csv')
van_df = van_df.drop(van_df.columns[0], axis=1)
van_df.Purchase = van_df.Purchase.map({'Yes':1, 'No': 0})
van_df.head()

Out[107]:

	MOSTYPE	MAANTHUI	MGEMOMV	MGEMLEEF	MOSHOOFD	MGODRK	MGODPR	MGODOV	MGODGE	MRELGE	...	ABRAND
0	33	1	3	2	8	0	5	1	3	7	...	1
1	37	1	2	2	8	1	4	1	4	6	...	1
2	37	1	2	2	8	0	4	2	4	3	...	1
3	9	1	3	3	3	2	3	2	4	5	...	1
4	40	1	4	2	10	1	4	1	4	7	...	1

5 rows × 86 columns

(a) Create a training set consisting of the first 1,000 observations, and a test set consisting of the remaining observations.

In [113]:

train,test = model_selection.train_test_split(van_df, test_size=van_df.shape[0] - 1000)

(b) Fit a boosting model to the training set with Purchase as the response and the other variables as predictors. Use 1,000 trees, and a shrinkage value of 0.01. Which predictors appear to be the most important?

In [146]:

train_x = train.drop('Purchase', axis=1)
train_y = train.Purchase
clr = ensemble.GradientBoostingClassifier(n_estimators=1000, learning_rate=0.01).fit(train_x, train_y)
test_x = test.drop('Purchase', axis=1)
test_y = test.Purchase
preds = clr.predict(test_x)
error = 1 - ((preds == test_y).sum() / preds.shape[0])

print('classification error: {}'.format(error))


# plot
_,_ = plt.subplots(figsize=(15, 7.5))
bar_df = pd.DataFrame({'predictor': train_x.columns, 'importance' : clr.feature_importances_})
sns.barplot(x='predictor', y='importance', data=bar_df.sort_values('importance', ascending=False).iloc[0:10])

classification error: 0.06574035669846534

Out[146]:

<matplotlib.axes._subplots.AxesSubplot at 0x11b7680f0>

(c) Use the boosting model to predict the response on the test data. Predict that a person will make a purchase if the estimated probability of purchase is greater than 20%. Form a confusion matrix. What fraction of the people predicted to make a purchase do in fact make one? How does this compare with the results obtained from applying KNN or logistic regression to this data set?

In [163]:

preds = (clr.predict_proba(test_x)[:,1] > 0.2)
conf_mtrx = pd.DataFrame(metrics.confusion_matrix(test_y, preds))

predicted_true = conf_mtrx[0][1] + conf_mtrx[1][1]
print('Positive Predictive Value: {}'.format(conf_mtrx[1][1] / predicted_true))

Positive Predictive Value: 0.1840277777777778

In [169]:

logit = linear_model.LogisticRegression().fit(train_x, train_y)
preds = (logit.predict_proba(test_x)[:,1] > 0.2)
conf_mtrx = pd.DataFrame(metrics.confusion_matrix(test_y, preds))
predicted_true = conf_mtrx[0][1] + conf_mtrx[1][1]
print('Positive Predictive Value: {}'.format(conf_mtrx[1][1] / predicted_true))

Positive Predictive Value: 0.1840277777777778

In [174]:

knn = neighbors.KNeighborsClassifier(n_neighbors=5).fit(train_x, train_y)
preds = (knn.predict_proba(test_x)[:,1] > 0.2)
conf_mtrx = pd.DataFrame(metrics.confusion_matrix(test_y, preds))
predicted_true = conf_mtrx[0][1] + conf_mtrx[1][1]
print('Positive Predictive Value: {}'.format(conf_mtrx[1][1] / predicted_true))

Positive Predictive Value: 0.0625

Apply boosting, bagging, and random forests to a data set of your choice. Be sure to fit the models on a training set and to evaluate their performance on a test set. How accurate are the results compared to simple methods like linear or logistic regression? Which of these approaches yields the best performance?

In [ ]:

# TODO: Come back for this one with a kaggle dataset and take XGBoost out for a spin

	MOSTYPE	MAANTHUI	MGEMOMV	MGEMLEEF	MOSHOOFD	MGODRK	MGODPR	MGODOV	MGODGE	MRELGE	...	ABRAND
0	33	1	3	2	8	0	5	1	3	7	...	1
1	37	1	2	2	8	1	4	1	4	6	...	1
2	37	1	2	2	8	0	4	2	4	3	...	1
3	9	1	3	3	3	2	3	2	4	5	...	1
4	40	1	4	2	10	1	4	1	4	7	...	1

	MOSTYPE	MAANTHUI	MGEMOMV	MGEMLEEF	MOSHOOFD	MGODRK	MGODPR	MGODOV	MGODGE	MRELGE	...	ABRAND
0	33	1	3	2	8	0	5	1	3	7	...	1
1	37	1	2	2	8	1	4	1	4	6	...	1
2	37	1	2	2	8	0	4	2	4	3	...	1
3	9	1	3	3	3	2	3	2	4	5	...	1
4	40	1	4	2	10	1	4	1	4	7	...	1

	MOSTYPE	MAANTHUI	MGEMOMV	MGEMLEEF	MOSHOOFD	MGODRK	MGODPR	MGODOV	MGODGE	MRELGE	...	ABRAND
0	33	1	3	2	8	0	5	1	3	7	...	1
1	37	1	2	2	8	1	4	1	4	6	...	1
2	37	1	2	2	8	0	4	2	4	3	...	1
3	9	1	3	3	3	2	3	2	4	5	...	1
4	40	1	4	2	10	1	4	1	4	7	...	1