CS-109A Introduction to Data Science

Lab 6: Logistic Regression¶

Harvard University
Fall 2019
Instructors: Pavlos Protopapas, Kevin Rader, Chris Tanner
Lab Instructors: Chris Tanner and Eleni Kaxiras.
Contributors: Will Claybaugh, David Sondak, Chris Tanner

Learning Goals¶

In this lab, we'll explore different models used to predict which of several labels applies to a new datapoint based on labels observed in the training data.

By the end of this lab, you should:

• Be familiar with the sklearn implementations of
  • Linear Regression
  • Logistic Regression
• Be able to make an informed choice of model based on the data at hand
• (Bonus) Structure your sklearn code into Pipelines to make building, fitting, and tracking your models easier
• (Bonus) Apply weights to each class in the model to achieve your desired tradeoffs between discovery and false alarm in various classes

Part 1: The Wine Dataset¶

The dataset contains 11 chemical features of various wines, along with experts' rating of that wine's quality. The quality scale technically runs from 1-10, but only 3-9 are actually used in the data.

Our goal will be to distinguish good wines from bad wines based on their chemical properties.

Read-in and checking¶

We do the usual read-in and verification of the data:

Building the training/test data¶

As usual, we split the data before we begin our analysis.

Today, we take the 'quality' variable as our target. There's a debate to be had about the best way to handle this variable. It has 10 categories (1-10), though only 3-9 are used. While the variable is definitely ordinal- we can put the categories in an order everyone agrees on- the variable probably isn't a simple numeric feature; it's not clear whether the gap between a 5 and a 6 wine is the same as the gap between an 8 and a 9.

Ordinal regression is one possibility for our analysis (beyond the scope of this course), but we'll view the quality variable as categorical. Further, we'll simplify it down to 'good' and 'bad' wines (quality at or above 7, and quality at or below 6, respectively). This binary column already exists in the data, under the name 'good'.

Now that we've split, let's explore some patterns in the data

It looks like there aren't any particularly strong correlations among the predictors (maybe sulfur dioxide and free sulfur dioxide) so we're safe to keep them all. It also looks like the different quality categories have roughly the same distribution of most variables, with volatile/fixed acidity and alcohol seeming like the most useful predictors.

Part 2 (Introduction): Binary Logistic Regression¶

Linear regression is usually a good baseline model, but since the outcome we're trying to predict only takes values 0 and 1 we'll want to use logistic regression instead of basic linear regression.

We'll begin with statsmodels, because cs109 likes confidence intervals and checking that coefficients make sense.

Let's talk about the output:
First, "optimization terminated successfully". Recall that linear regression and its simple formula for the optimal betas is a rarity in machine learning and statistics: most models are fit to the data algorithmically, not via a formula. This message is letting us know that the algorithm seems to have worked.

Second, the pseudo $R^2$ is rather low (.23). As with regular $R^2$, we might take this as a sign that the model is struggling.

Finally, let's look at the coefficients.

• Several of the coefficients are statistically significant, including
  • Fixed acidity - good
  • Volatile Acidity - bad
  • Residual Sugar - good (judge have a sweet tooth?)
  • Chlorides - bad
  • Sulphates - good
  • Alcohol - good (judges like getting drunk?)
  • The rest only reach a coefficient size we would often observe by chance alone, without any actual effect from the predictor

More formal interpretations are of coefficients are long-winded. "A one unit increase in alcohol (holding all else constant) results in a predicted 0.494 increase in the log odds of a wine being classified as good".

We can't be more precise because the effect of one unit of alcohol depends on how much alcohol there already is. The one unit increase/decrease matters more if the wine is otherwise on the border between good and bad. If the wine is undrinkable (in the far left tail of the sigmoidal curve) one unit of alcohol barely moves the probability, while if the wine is in the middle of the curve that unit of acidity has much more practical impact.

1. Are there any bones you'd like to pick with the model I've laid out? Can you think of a better logistic regression model?

Prediction¶

One of the really cool features of logistic regression is that it hands back probabilities of a given case being 1 or 0, rather than just 1s and 0s. That lets us do neat things like set different cutoffs for what counts as a 1 and do ROC analysis and so on. Here, we'll just set the cutoff at 0.5: if a 1 is reported as more likely, predict a 1. (You can play with the cutoff yourself and see if you can make the model do better by trading false positives and false negatives).

Because this is statsmodels, we'll need to import a tool or do the test set score calculation ourselves. Here, it's easy enough to implement:

• do the predictions
• compare with .5
• find out what percentage of our binary predictions matched the truth

Wow! 80% is a pretty good performance! We can pretty much tell the bad wines from the good.

Here's a sanity check:

Oh... no... wait. A model that says "all wines are bad" also scores 80% on the test set. Our fancy model isn't really doing that well.

Moral of the story: Before you congratulate a model, think of a truly trivial model to compare it to.

1. Re-create the results above but this time work with sklearn. Use the LogisticRegression class. Follow the usual .fit, .score procedure. To match statsmodel's coefficient values (roughly), you will need to adjust the input parameters:
   • C
   • solver
   • One other parameter
   • See the sklearn documentation

Hint: statsmodels uses a Newton-Raphson method to optimize the beta values.

Answer:

Speaker note: When presenting solution, model reading the documentation from the webpage. How does one know where to look?
Speaker note: Mention the wide variety of solvers and how (some) use different levels of derivatives to converge in fewer steps

The Decision Boundary¶

One powerful way to think about classification models is to consider where and how they draw the line between predicting "class A" and "class B". The code below lets you play with a 2d logistic regression. Points towards yellow will be predicted as 1s, points towards violet will be predicted as 0s.

1. What happens to the decision boundary as the coefficient on X increases?
2. What happens if you increase the Y coefficient to match?
3. What does the constant term control?
4. What impact do higher-order and interaction terms have on the boundary?
5. What parameter settings should I show the class?

Answers:

your answer here

1. The boundary tips towards vertical
2. The boundary is in the same place as it was originally, but is squished together. The model is much more certain about how to predict points a given distance from the boundary
3. It shifts the boundary, perpendicular to its current orientation
4. Including squared terms allows quadratic decision boundaries, and the interraction term allows hyperbolic boundaries

Part 3 (The Real Challenge): Multiclass Classification¶

Before we move on, let's consider a more common use case of logistic regression: predicting not just a binary variable, but what level a categorical variable will take. Instead of breaking the quality variable into 'good' and 'other', let's discretize into 'good, 'medium', and 'bad'.

The cut function obviously stores a lot of extra information for us. It's a very useful tool for discretizing an existing variable.

1. Adapt your earlier logistic regression code to fit to the new training data. What is stored in .coef_ and .intercept_?
2. How well does this model predict the test data?
3. Put the model's performance in context. Think of a trivial model to compare to, and provide its accuracy score on the test set.

Answers:

1.

your answer here

1. We get three sets of coefficients, and three intercepts.

We need three sets because (under the default 'one versus rest' strategy) we fit three models. When predicting, model 1 reports a probability of the new example coming from class A or from the cloud of remaining classes. Model 2 reports the probability of whether the example comes from class B or the cloud of remaining classes, and so on. We take this set of scores and pick the biggest one (we classify as whichever class has the biggest ratio of "this class" to "not this class").

2.

your answer here

2. The model does pretty well at predicting the test data...

3.

your solution here

But, a trivial model that guesses the most likely class also does really well on the test set, too.

Summary¶

• Logistic regression extends OLS to work naturally with a dependent variable that's only ever 0 and 1.
• In fact, Logistic regression is even more general and is used for predicting the probability of an example belonging to each of $N$ classes.
• The code for the two cases is identical and just like LinearRegression: .fit, .score, and all the rest
• Significant predictors does not imply that the model actually works well. Signifigance can be driven by data size alone.
• The data aren't enough to do what we want

Warning: Logistic regression tries to hand back valid probabilities. As with all models, you can't trust the results until you validate them- if you're going to use raw probabilities instead of just predicted class, take the time to verify that if you pool all cases where the model says "I'm 30% confident it's class A" 30% of them actually are class A.

Part 4: Dimensionality Reduction¶

Our models are clearly struggling, but it's hard to tell why. Let's PCA to shrink the problem down to 2d (with as little loss as possible) and see if that gives us a clue about what makes this problem tough.

Some comments:

1. Both scaling and reducing dimension follow the same pattern: we fit the object to the training data, then use .transform to convert the training and test data. This ensures that, for instance, we scale the test data using the training mean and variance, not its own mean and variance
2. We need to equalize the variance of each feature before applying PCA, otherwise certain dimensions will dominate the scaling: our PCA dimensions would just be the features with the largest spread.

Well, that gives us some idea of why the problem is difficult: the good wines and bad wines are hiding right among the average wines. It does look like the wines separate into two groups, though, possibly one for reds and one for whites.

1. What critique can you make against the plot above? Why does this plot not prove that the different wines are hopelessly similar?
2. The wine data we've used so far consist entirely of continuous predictors. Would PCA work with categorical data?

Answer:

your answer here

1. The PCA dimensions are chosen without regard to the y variable. Thus it is possible that the very next PCA dimension will lift the red points up out of the page, push the blue points down into it, and leave the cyan points where they are; such a dimension would separate the different types of wine and make classification easy.
2. PCA would not work with categorical data. PCA requires there to be a meaningful notion of distance between points. Categorical or ordinal data is not enough.

1. Edit the code above to plot the locations of red wines and white wines

Answer:

Evaluating PCA - Variance Explained¶

One of the criticisms we made of the PCA plot was that it's lost something from the original data.

Let's actually investigate how much of the original data's structure the 2d PCA captures. We'll look at the explained_variance_ratio_ portion of the PCA fit. This lists, in order, the percentage of the x data's total variance that is captured by the nth PCA dimension.

The first PCA dimension captures 33% of the variance in the data, and the second PCA dimension adds another 20%. Together, we've got about half of the total variation in the training data covered with just these two dimensions.

1. Fit a PCA that finds the first 10 PCA components of our training data
2. Use np.cumsum to print out the variance we'd be able to explain by using n PCA dimensions for n=1 through 10
3. Does the 10-dimension PCA agree with the 2d PCA on how much variance the first components explain? Do the 10d and 2d PCAs find the same first two dimensions? Why or why not?
4. Make a plot of number of PCA dimensions against total variance explained. What PCA dimension looks good to you?

Hint: np.cumsum stands for 'cumulative sum', so np.cumsum([1,3,2,-1,2]) is [1,4,6,5,7]

Answer:

3.

your answer here

The 10d PCA and the 2d PCA agree about how much variance the first two components explain. The 10d and 2d PCA give the same components in the same order. This means it's safe to simply fit a PCA with the largest number of components you expect you will need, and take a subset as appropriate.

4.

A PCA dimension of 3, 4, or 5 looks good to me. These values are roughly where we hit diminishing returns on variance explained.

Plots like the one above are called 'Scree' or 'Elbow' plots. They are often used to heuristically select a good number of PCA dimensions.

Summary¶

• PCA maps a high-dimensional space into a lower dimensional space.
• The PCA dimensions are ordered by how much of the original data's variance they capture
  • There are other cool and useful properties of the PCA dimensions (orthogonal, etc.). See a textbook.
• PCA on a given dataset always gives the same dimensions in the same order.
• You can select the number of dimensions by fitting a big PCA and examining a plot of the cumulative variance explained.

Part 5: Did we fail?¶

None of the models worked, and we can't tell good wines from bad. Was it all a waste of time and money?

Not really. All analyses are a roll of the dice. Some analyses fail, like this one did, becuase the data at hand just don't support the task we've set out.

What can we do about it?

1. Be honest about the methods and the null result. Lots of analyses fail.
2. Collect a dataset you think has a better chance of success. Maybe we collected the wrong chemical signals...
3. Keep trying new approaches. Just beware of overfitting the data you're validating on. Always have a test set locked away for when the final model is built.
4. Change the question. Maybe something you noticed during analysis seems interesting or useful (classifying red versus white). But again, you the more you try, the more you might overfit, so have test data locked away.
5. Just move on. If the odds of success start to seem small, maybe you need a new project.

The Moral of the Lab¶

• Sometimes, the data just aren't enough to adequately predict outcomes.
• In this lab we saw that no amount of modeling finesse would let us use a wine's chemical properties to tell good wines and bad wines from mediocre ones.
  • The chemical properties were very good at telling red wines from whites, however.
• PCA helped us visualize the data and confirm that the highly rated wines just aren't chemically distinct from the other wines.
• NOT ALL ANALYSES YIELD USEFUL RESULTS Sometimes (arguably most of the time), the data aren't suitable for a task or just don't have anything interesting to say.

Part 6 (Sidebar): Pipelines¶

Remember when we were trying to adapt our LDA model to run on the training data with 'red' dropped? We had to invent new variable names and define functions and it was generally much harder than it needed to be. Pipelines are sklearn's tool for packaging an entire analysis together into a single object. This enables convenient inspection, saving, deployment, and (yes) cross validation of the model.

Let's look at an example (we'll switch the model to KNN to justify some later analysis).

There's also a convenience function make_pipeline that lets us skip naming the different steps. Notice the default names are all-lowercase versions of the class names (standardscaler, pca, kneighborsclassifier)

It's easy to run the whole modelling process on new data:

As promised, cross validation tools work directly with the pipeline object.

Note: In general, more PCA components will work better for prediction. However, KNN often performs worse as dimension increases, meaning there may be a meaningful balance point between capturing more variance and a space small enough for KNN to work well.

Part 7 (Sidebar): Weighting the training points¶

Some models can accept weights on the training points to given them greater priority in the model's fitting process. This can be useful if, for instance, certain classes are rare but we want to be sure the model classifies them correctly (e.g. we're trying to classify cancers and one form is rare but very aggressive). In general, weighting training points is like moving along the ROC curve; we change some model parameters to alter the mistakes the model makes to be more in line with our tastes.

Let's see this in action with a logistic regression:

Without weighting, the model plays it safe and predicts that all of the test set wines are medium. With weighting, the model is told to care more about getting the bad and good wines right. The model does as we've asked and correctly IDs 3 good/bad test wines, at the price of 17 falsely bad wines and 16 falsely good wines. However, if identifying bad and good wines is, as implied, 100 times more important than identifying medium wines, we've made a really good trade.

1. What happens if you give a weight of 0 to the medium wines?
2. What weighting gives results that accord with your personal sense of what the model should be doing? How many actually-medium bottles is a single good bottle worth?

Answers:

1. The model learns a classification rule that never predicts 'medium'. It's as it we dropped the medium wines from training.
2. 100, 1, 100 looks the best to me. We get a 1-in-8 sucess rate on the wines flagged as good. However, I found these values by looking at the test set confusion matrix; it's not clear they'd maintain the 1-in-8 ratio on new data.