Classification

weekly = sm.datasets.get_rdataset("Weekly", package="ISLR")
print(weekly.data.info())

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 1089 entries, 0 to 1088
Data columns (total 9 columns):
 #   Column     Non-Null Count  Dtype
---  ------     --------------  -----
 0   Year       1089 non-null   int64
 1   Lag1       1089 non-null   float64
 2   Lag2       1089 non-null   float64
 3   Lag3       1089 non-null   float64
 4   Lag4       1089 non-null   float64
 5   Lag5       1089 non-null   float64
 6   Volume     1089 non-null   float64
 7   Today      1089 non-null   float64
 8   Direction  1089 non-null   object
dtypes: float64(7), int64(1), object(1)
memory usage: 76.7+ KB
None

There is a non-numeric column Direction. We need some more detail about the data set to see if we can convert the non-numeric column to a numeric column.

print(weekly.data.describe())

              Year         Lag1         Lag2         Lag3         Lag4  \
count  1089.000000  1089.000000  1089.000000  1089.000000  1089.000000
mean   2000.048669     0.150585     0.151079     0.147205     0.145818
std       6.033182     2.357013     2.357254     2.360502     2.360279
min    1990.000000   -18.195000   -18.195000   -18.195000   -18.195000
25%    1995.000000    -1.154000    -1.154000    -1.158000    -1.158000
50%    2000.000000     0.241000     0.241000     0.241000     0.238000
75%    2005.000000     1.405000     1.409000     1.409000     1.409000
max    2010.000000    12.026000    12.026000    12.026000    12.026000

              Lag5       Volume        Today
count  1089.000000  1089.000000  1089.000000
mean      0.139893     1.574618     0.149899
std       2.361285     1.686636     2.356927
min     -18.195000     0.087465   -18.195000
25%      -1.166000     0.332022    -1.154000
50%       0.234000     1.002680     0.241000
75%       1.405000     2.053727     1.405000
max      12.026000     9.328214    12.026000

weekly.data["Direction"].unique()

array(['Down', 'Up'], dtype=object)

So we see that the Directions column has only two unique values, Up, and Down, showing if the market was up or down. If required we can easily convert this to a numeric column with the map Up: 1, Down: 0. For a graphical summary of the data, we use either pair plots or a heat map of the correlation matrix.

spm = sns.pairplot(weekly.data)
spm.fig.set_size_inches(12, 12)
spm.savefig("img/4.10.a_pair.png")

I am not really seeing any pattern in the pair plots except that the lags are always between -10 and 10, and there does not seem to be a lot correlation between the different lags. The correlation matrix will help to make this more concrete.

corr_mat = weekly.data[weekly.data.columns[:-1]].corr()

fig, ax = plt.subplots(figsize=(8, 6))
cmap = sns.diverging_palette(220, 10, sep=80, n=7)
mask = np.triu(np.ones_like(corr_mat, dtype=np.bool))

with sns.axes_style("white"):
    sns.heatmap(corr_mat, cmap=cmap, mask=mask, robust=True, annot=True, ax=ax)

plt.tight_layout()
fig.savefig("img/4.10.a_corr_mat.png", dpi=90)
plt.close()

This matches the observation from the pair plots. There is a high correlation between Year and Volume, but everything else is mostly uncorrelated with each other.

(Note: Constructing seaborn pair plots are somewhat computationally intensive. Might be better to start with the correlation matrix.)

We will train a logistic regression model on the entire data set with Direction as response and Lag1-5, and Volume as predictors. We will use the logistic regression implementation from statsmodels. There are two ways of doing it, using statsmodels.api or statsmodels.formula.api. I like the elegant R-like formula syntax that statsmodels.formula.api provides, and I am not a big fan of the exog / endog nomenclature used by statsmodels.api (endog, exog, what’s that?). I will be using the statsmodels.formula.api wherever possible.

More importantly we need to encode Direction as discussed earlier for statsmodels logistic regression, scikit-learn can work with categorical variables.

weekly_data = weekly.data.copy()
weekly_data["Direction"] = weekly.data["Direction"].map({"Up": 1, "Down": 0})

logit_model = smf.logit("Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 + Volume",
                        data=weekly_data).fit()
print(logit_model.summary())

Optimization terminated successfully.
         Current function value: 0.682441
         Iterations 4
                           Logit Regression Results
==============================================================================
Dep. Variable:              Direction   No. Observations:                 1089
Model:                          Logit   Df Residuals:                     1082
Method:                           MLE   Df Model:                            6
Date:                Thu, 18 Jun 2020   Pseudo R-squ.:                0.006580
Time:                        18:11:50   Log-Likelihood:                -743.18
converged:                       True   LL-Null:                       -748.10
Covariance Type:            nonrobust   LLR p-value:                    0.1313
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept      0.2669      0.086      3.106      0.002       0.098       0.435
Lag1          -0.0413      0.026     -1.563      0.118      -0.093       0.010
Lag2           0.0584      0.027      2.175      0.030       0.006       0.111
Lag3          -0.0161      0.027     -0.602      0.547      -0.068       0.036
Lag4          -0.0278      0.026     -1.050      0.294      -0.080       0.024
Lag5          -0.0145      0.026     -0.549      0.583      -0.066       0.037
Volume        -0.0227      0.037     -0.616      0.538      -0.095       0.050
==============================================================================

Among the predictor variable only the Lag2 appear to be statistically significant. Even then the p-value of 0.030 it still relatively large as compared to p-values of statistically significant predictors that we have seen in previous chapters. So Lag2 might still turn out to be statistically insignificant.

I got the idea of using percentages and labels in the heat map of the confusion matrix from Confusion Matrix Visualization.

def make_confusion_matrix_heatmap(conf_mat, categories, ax):
    """
    Makes a heat map visualization of the confusion matrix.
    """
    group_names = ["True Neg", "False Pos", "False Neg", "True Pos"]
    group_counts = [f"{value:.0f}" for value in conf_mat.flatten()]
    group_percentages = [f"{value:.2%}" for value in
                         conf_mat.flatten()/np.sum(conf_mat)]
    labels = [f"{v1}\n{v2}\n{v3}" for v1, v2, v3 in
              zip(group_names,group_counts,group_percentages)]
    labels = np.asarray(labels).reshape(2,2)

    with sns.axes_style("white"):
        sns.heatmap(conf_mat, cmap="Blues", fmt="", annot=labels, cbar=False,
                    xticklabels=categories, yticklabels=categories, ax=ax)

conf_mat = logit_model.pred_table(0.5)

fig, ax = plt.subplots(figsize=(4, 4))
categories = ["Down", "Up"]
make_confusion_matrix_heatmap(conf_mat, categories=categories, ax=ax)
plt.tight_layout()
fig.savefig("img/4.10.c_conf_mat.png", dpi=90)
plt.close()

The values along the diagonal give the number of correct predictions. There were 54 instances of both the predicted and the observed responses being Down, and 557 instances of both the predicted and the observed responses being Up . There were 430 instances where the observed response was Down, but the model predicted it to be Up. These are the false positives. And there were 48 instances where the observed response was Up, but the model predicted it to be Down. The accuracy of the model is defined as the ratio of the total number of correct predictions to the total number of predictions.

accuracy = (np.sum(np.diag(conf_mat))) / np.sum(conf_mat)
print(f"Accuracy: {accuracy:.2%}")

Accuracy: 56.11%

This is essentially the sum of the percentages on the diagonals in the above confusion matrix heat map.

Now we will train the logistic regression model to predict Direction with only the data from 1990 to 2008, and with Lag2 as the only predictor.

weekly_training_set = weekly_data[weekly_data["Year"].between(1990, 2008)]
weekly_test_set = weekly_data.drop(weekly_training_set.index)

logit_model2 = smf.logit("Direction ~ Lag2", data=weekly_training_set).fit()
print(logit_model2.summary())

Optimization terminated successfully.
         Current function value: 0.685555
         Iterations 4
                           Logit Regression Results
==============================================================================
Dep. Variable:              Direction   No. Observations:                  985
Model:                          Logit   Df Residuals:                      983
Method:                           MLE   Df Model:                            1
Date:                Thu, 18 Jun 2020   Pseudo R-squ.:                0.003076
Time:                        19:22:59   Log-Likelihood:                -675.27
converged:                       True   LL-Null:                       -677.35
Covariance Type:            nonrobust   LLR p-value:                   0.04123
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept      0.2033      0.064      3.162      0.002       0.077       0.329
Lag2           0.0581      0.029      2.024      0.043       0.002       0.114
==============================================================================

Here we want the confusion matrix for the test set. The builtin pred_table method from statsmodels does not work with data that has not been used for training. We will use the confusion_matrix method from sklearn.metrics. An alternative is to use numpy directly, as shown here.

pred = np.array(logit_model2.predict(weekly_test_set["Lag2"]) > 0.5, dtype="int")
conf_mat2 = confusion_matrix(weekly_test_set["Direction"], pred)

fig, ax = plt.subplots(figsize=(4, 4))
categories = ["Down", "Up"]
make_confusion_matrix_heatmap(conf_mat2, categories=categories, ax=ax)
plt.tight_layout()
fig.savefig("img/4.10.d_conf_mat.png", dpi=90)
plt.close()

We can see that the accuracy of this model is more (62.50%), even when it has not been trained on the test data set.

We use the same training and test data sets as before but now to fit a Linear discriminant analysis model. We will use the linear discriminant analysis implementation from scikit-learn.

lda_model = LinearDiscriminantAnalysis()
lda_model.fit(weekly_training_set["Lag2"].values[:, None],
              weekly_training_set["Direction"].values)

preds = lda_model.predict(weekly_test_set["Lag2"].values[:, None])
conf_mat_lda = confusion_matrix(weekly_test_set["Direction"], preds)

fig, ax = plt.subplots(figsize=(4, 4))
categories =["Down", "Up"]
make_confusion_matrix_heatmap(conf_mat_lda, categories, ax)
plt.tight_layout()
fig.savefig("img/4.10.e_conf_mat.png", dpi=90)
plt.close()

The LDA model has practically the same accuracy as the previous logistic regression model.

We will repeat the above for the quadratic discriminant analysis model.

qda_model = QuadraticDiscriminantAnalysis()
qda_model.fit(weekly_training_set["Lag2"].values[:, None],
              weekly_training_set["Direction"].values)

preds = qda_model.predict(weekly_test_set["Lag2"].values[:, None])
conf_mat_qda = confusion_matrix(weekly_test_set["Direction"], preds)

fig, ax = plt.subplots(figsize=(4, 4))
categories =["Down", "Up"]
make_confusion_matrix_heatmap(conf_mat_qda, categories, ax)
plt.tight_layout()
fig.savefig("img/4.10.f_conf_mat.png", dpi=90)
plt.close()

The QDA model does not give any negative results, but now the percentage of false positive rivals the percentage of true positives. The accuracy has decreased to 58.65%.

Rinse and repeat but now with KNN for \(K = 1\).

knn_model = KNeighborsClassifier(n_neighbors=1)
knn_model.fit(weekly_training_set["Lag2"].values[:, None],
              weekly_training_set["Direction"].values)

preds = knn_model.predict(weekly_test_set["Lag2"].values[:, None])
conf_mat_knn = confusion_matrix(weekly_test_set["Direction"], preds)

fig, ax = plt.subplots(figsize=(4, 4))
categories =["Down", "Up"]
make_confusion_matrix_heatmap(conf_mat_knn, categories, ax)
plt.tight_layout()
fig.savefig("img/4.10.g_conf_mat.png", dpi=90)
plt.close()

The KNN model with \(K = 1\) has an accuracy of 49.04%.

Logistic regression and Linear discriminant analysis appear to provide the best results on this data as they have the highest accuracy among the four models. KNN with \(K = 1\) gave the worst performance.

auto_data = sm.datasets.get_rdataset("Auto", "ISLR").data

median_mileage = auto_data["mpg"].median()
auto_data["mpg01"] = np.where(auto_data["mpg"] > median_mileage, 1, 0)
print(tabulate(auto_data.head(), auto_data.columns, tablefmt="orgtbl"))

|    |   mpg |   cylinders |   displacement |   horsepower |   weight |   acceleration |   year |   origin | name                      |   mpg01 |
|----+-------+-------------+----------------+--------------+----------+----------------+--------+----------+---------------------------+---------|
|  1 |    18 |           8 |            307 |          130 |     3504 |           12   |     70 |        1 | chevrolet chevelle malibu |       0 |
|  2 |    15 |           8 |            350 |          165 |     3693 |           11.5 |     70 |        1 | buick skylark 320         |       0 |
|  3 |    18 |           8 |            318 |          150 |     3436 |           11   |     70 |        1 | plymouth satellite        |       0 |
|  4 |    16 |           8 |            304 |          150 |     3433 |           12   |     70 |        1 | amc rebel sst             |       0 |
|  5 |    17 |           8 |            302 |          140 |     3449 |           10.5 |     70 |        1 | ford torino               |       0 |

We want to graphically see how does mpg01 associate with the other features in our modified Auto data set. For this purpose we will use correlation matrix and boxplots. We can also use paired scatterplots (pairplots in seaborn parlance), but as I mentioned earlier seaborn's implementation of pairplots is quite computationally intensive and so far I have not seen it provide any information for classification problems that correlation matrix does not.

corr_mat = auto_data[auto_data.columns[1:]].corr()

cmap = "RdBu"
mask = np.triu(np.ones_like(corr_mat, dtype=np.bool))

fig, ax = plt.subplots(figsize=(8, 6))
with sns.axes_style("white"):
    sns.heatmap(corr_mat, cmap=cmap, mask=mask, annot=True, robust=True, ax=ax)

plt.tight_layout()
fig.savefig("img/4.11.b_corr_mat.png", dpi=90)
plt.close()

We see that mpg01 has high negative correlation with cylinders, displacement, horsepower, and weight. These features also high correlation with each other, for example displacement and cylinders have a correlation of 0.95. This means if our model includes displacement then we do not gain any new information by adding cylinders to it. Presumably we will be able to predict mpg01 with just horsepower, one of displacement, cylinders and weight, and maybe origin.

Let us now see what boxplots tell us.

fig, axs = plt.subplots(4, 2, figsize=(12, 10))
sns.boxplot(y="cylinders", x="mpg01", data=auto_data, ax=axs[0, 0])
sns.boxplot(y="displacement", x="mpg01", data=auto_data, ax=axs[0, 1])
sns.boxplot(y="horsepower", x="mpg01", data=auto_data, ax=axs[1, 0])
sns.boxplot(y="weight", x="mpg01", data=auto_data, ax=axs[1, 1])
sns.boxplot(y="acceleration", x="mpg01", data=auto_data, ax=axs[2, 0])
sns.boxplot(y="year", x="mpg01", data=auto_data, ax=axs[2, 1])
sns.boxplot(y="origin", x="mpg01", data=auto_data, ax=axs[3, 0])
plt.tight_layout()
fig.savefig("img/4.11.b_box_plots.png")
plt.close()

We see that having more cylinders (or displacement or horsepower or weight) generally gives you low mileage. The origin really does not matter, since except for a few outliers all have a mpg01 of 1. And as far as year and acceleration are concerned, there are a fair bit of overlap between the mpg01 = 0 and mpg01 = 1 groups, and it is difficult to draw any clear conclusion from them. This essentially corroborates what we had learned from the correlation matrix.

We will use scikit-learn's train_test_split to create training and test sets.

X = auto_data[["cylinders", "displacement", "horsepower", "weight"]]
y = auto_data["mpg01"]
train_X, test_X, train_y, test_y = train_test_split(X, y, random_state=42)

The question asks us to use all the features that seemed most associated with mpg01. Even though I mentioned that the high correlation between cylinders, displacement, and weight means that only one of them will be sufficient, we will still use all of them here.

lda_model = LinearDiscriminantAnalysis()
lda_model.fit(train_X, train_y)

test_pred = lda_model.predict(test_X)
conf_mat = confusion_matrix(test_pred, test_y)

fig, ax = plt.subplots(figsize=(4, 4))
make_confusion_matrix_heatmap(conf_mat, [0, 1], ax)
fig.savefig("img/4.11.d_conf_mat.png", dpi=90)
plt.close()

This model has a very high accuracy of 88.78%. Or in other words the test error is ≈ 0.11. Just for curiosity I want to try a model with just cylinders as the predictor.

train_X_cyl = train_X["cylinders"]
test_X_cyl = test_X["cylinders"]

lda_model_cyl = LinearDiscriminantAnalysis()
lda_model_cyl.fit(train_X_cyl.values[:, None], train_y)

test_pred_cyl = lda_model_cyl.predict(test_X_cyl.values[:, None])
conf_mat_cyl = confusion_matrix(test_pred_cyl, test_y)

fig, ax = plt.subplots(figsize=(4, 4))
make_confusion_matrix_heatmap(conf_mat_cyl, [0, 1], ax)
fig.savefig("img/4.11.d_conf_mat_cyl.png", dpi=90)
plt.close()

This model has a comparable accuracy of 87.76% (test error ≈ 0.12), the true positives decreased by 1, while the false positives increased by 1. Essentially we got the same result with a single feature as we did by using four features. So the other features are indeed superfluous.

Following the realization made in the previous question, I will use only cylinders for prediction using quadratic discriminant analysis.

qda_model_cyl = QuadraticDiscriminantAnalysis()
qda_model_cyl.fit(train_X_cyl.values[:, None], train_y)

test_pred_cyl = qda_model_cyl.predict(test_X_cyl.values[:, None])
conf_mat_cyl_qda = confusion_matrix(test_pred_cyl, test_y)

fig, ax = plt.subplots(figsize=(4, 4))
make_confusion_matrix_heatmap(conf_mat_cyl_qda, [0, 1], ax)
fig.savefig("img/4.11.d_conf_mat_cyl_qda.png", dpi=90)
plt.close()

This model has the same accuracy as the LDA model with just cylinders. For completeness I will also try a second model with the other features.

qda_model = QuadraticDiscriminantAnalysis()
qda_model.fit(train_X, train_y)

test_pred = qda_model.predict(test_X)
conf_mat_qda = confusion_matrix(test_pred, test_y)

fig, ax = plt.subplots(figsize=(4, 4))
make_confusion_matrix_heatmap(conf_mat_qda, [0, 1], ax)
fig.savefig("img/4.11.d_conf_mat_qda.png", dpi=90)
plt.close()

The accuracy did not change.

Since we are not interested in the summary of the logistic regression model, we will use scikit-learn's implementation of logistic regression. In fact this will be template by which we decide whether to use statsmodels or scikit-learn: If we are more interested in the statistical summary of the model then we will statsmodels. If we are more interested in predictions then we will use scikit-learn.

logit_model = LogisticRegression()
logit_model.fit(train_X_cyl.values[:, None], train_y)

test_pred = logit_model.predict(test_X_cyl.values[:, None])
conf_mat_logit = confusion_matrix(test_pred, test_y)

fig, ax = plt.subplots(figsize=(4, 4))
make_confusion_matrix_heatmap(conf_mat_logit, [0, 1], ax)
fig.savefig("img/4.11.d_conf_mat_logit.png", dpi=90)
plt.close()

Same accuracy as the QDA model.

We will try KNN with \(K ∈ \{1, 10, 100\}\).

for k in np.logspace(0, 2, num=3, dtype=int):
    knn_model = KNeighborsClassifier(n_neighbors=k)
    knn_model.fit(train_X_cyl.values[:, None], train_y)
    acc = accuracy_score(knn_model.predict(test_X_cyl.values[:, None]), test_y)
    print(f"K: {k:3}, Accuracy: {acc:.2%}, Test error: {1 - acc:.2f}")

K:   1, Accuracy: 87.76%, Test error: 0.12
K:  10, Accuracy: 87.76%, Test error: 0.12
K: 100, Accuracy: 87.76%, Test error: 0.12

We get the same accuracy for all the three values of \(K\). We will try this again including the other features.

for k in np.logspace(0, 2, num=3, dtype=int):
    knn_model = KNeighborsClassifier(n_neighbors=k)
    knn_model.fit(train_X, train_y)
    acc = accuracy_score(knn_model.predict(test_X), test_y)
    print(f"K: {k:3}, Accuracy: {acc:.2%}, Test error: {1 - acc:.2f}")

K:   1, Accuracy: 83.67%, Test error: 0.16
K:  10, Accuracy: 85.71%, Test error: 0.14
K: 100, Accuracy: 85.71%, Test error: 0.14

We see that for \(K = 10, 100\) the accuracy is better compared to the \(K = 1\) case. However the accuracy has in fact decreased from the previous model with just cylinder as the predictor. This is most likely due to the "curse of dimensionality". This an empirical "proof" that for KNN low-dimensional models, whenever possible, are better.