Remember how we noted a couple of things down on our to do list to check for later? We said that we'd check 1
for a potential problem in our regression. 2
And the problem that we were concerned about was called multicollinearity. 3
The reason was that we had some high correlations between our features, 4
so multicollinearity was something that we wanted to test formally. 5
Now put simply multicollinearity is when two or more predictor variables in a regression are highly 6
related to one another. 7
In that case they do not provide unique or independent information to the model and the consequences 8
of multicollinearity are as follows. 9
There is a loss of reliability in the estimates of the effects for the individual features, 10
in particular the features that are affected. 11
In other words, there is a high variability in the coefficient estimates for small changes in the model, 12
for example, adding or removing a feature can have dramatic effects. 13
The estimates of our theta values in our model become unstable and can even switch signs from positive 14
to negative and vice versa. 15
The third effect of the problem is that the findings are strange or misleading or don't make sense. 16
These are the main symptoms of the multicollinearity problem. 17
Now, given these symptoms we can definitely say that on the third one our regression definitely passed, 18
because we did a sense check on the signs of our coefficients and we found that the coefficients definitely 19
had the right signs. 20
They made sense from a logical perspective, 21
like more rooms increased the property prices and more pollution decreased the property prices for example. 22
But in this lesson I'm going to show you how you can do a formal check, how you can look at a metric to 23
tell you whether you have this problem or not. 24
The formal check that we will do is by looking at a statistic called the variance inflation factor or 25
VIF and the second thing that we will do is we will monitor if moving a feature causes a dramatic 26
change in our theta parameters, but that will wait until the next lesson. 27
Now let's go back to Jupyter notebook and add a new section heading. This section heading is gonna be 28
called "Testing for Multicollinearity", so "Testing for Multicollinearity", and to test for multicollinearity, 29
we said we're gonna be using a statistic called the variance inflation factor. 30
So let me quickly explain how the variance inflation factor does its job, because this statistic is a 31
measure of collinearity among the features within a multiple regression and what it will do is it 32
will spit out a number that quantifies the severity of multicollinearity, so in effect we get a number 33
similar to how we get a number with the p-values and we can look at this number and compare it to a 34
threshold, is it above the threshold or is it below the threshold. 35
That's how it's gonna work. 36
But before we write the Python code to actually calculate the variance inflation factors, let's add some 37
formulas in LaTeX notation to the section heading to kind of better understand what's going on. 38
I'm going to show you how the variance inflation factor is calculated for one of the features 39
as an example. The first thing that will happen is that a regression will be run on that feature. 40
Let's say we're looking at our TAX features. 41
So what will happen is that we will look to explain the TAX values as a linear combination of all the 42
other features in the data set. 43
So you're gonna get a regression that reads something like "TAX = " and then some intercept, say "\alpha 44
_0 ", plus some parameter times the number of rooms, 45
so " \alpha _1 RM", plus some other parameter "\alpha _2" times the next feature, 46
say nitrous oxide (NOX), plus and so on till you get to LSTAT, which will be "\alpha _ 47
{12} 48
LSTAT". Since we have 13 features total and we're looking to explain TAX on property prices, 49
we're left with 12 different kind of parameters in front of the other features. 50
I've also used alphas here in my LaTeX notation so we don't get confused with the theta parameters that 51
we used earlier. 52
Now this is only step one of calculating the variance inflation factor. 53
Step 2 looks like this. In step 2 we get the variance inflation factor for our TAX feature. 54
So "VIF _{TAX} = 55
\frac {1}/{(1-R _ 56
{TAX} ^ 2)}". Now let me quickly 57
delete this backslash here, we don't need that, and hit Shift+Enter to make this a little more clear. This 58
is step 2. In Step 1 a regression is being run of all the other features against TAX and in step 2 the 59
r-squared of that regression is used to calculate the variance inflation factor, so 1 divided by 60
(1 - r^2) of this regression up top is the variance inflation factor for the TAX feature, so 61
that r-squared here is the r-squared from trying to explain the tax feature in terms of all the other 62
features in the dataset. So the point I'm trying to make is that we can calculate a variance inflation 63
factor for all the different features, right, just by swapping out this TAX for RM we get the variance 64
inflation factor for RM, so in our Python code what we're gonna do is we're going to calculate 13 different 65
variants inflation factors. Now the statsmodels module actually makes this very, very easy. We're not 66
gonna have to run this like two step calculation manually, we simply call a function, but I wanted to 67
show you what's going on behind the scenes in the LaTeX notation nonetheless. 68
So let's scroll to the top of our notebook and actually import this functionality. I'm going to come here 69
and I'm going to write "from statsmodels.stats.outliers_influence import variance_ 70
inflation_factor". 71
So, I'm importing this specific function from "statsmodels.stats.outliers_influence" 72
Now let me hit Shift+Enter on the cell. 73
Interestingly, the deprecation warning disappeared when I did that. And I'm going to scroll back down and 74
now I'm actually going to call this function. 75
Check it out. 76
"variance_inflation_factor()" 77
and before I put any arguments between these two parentheses I'm going to hit Shift+Tab on my keyboard to 78
bring up the quick documentation on this function. 79
What we see here is that we need two inputs, exog and exog_idx, hit this little plus sign 80
here, 81
scroll down a little bit and we can see that the first argument is an n-dimensional array 82
and the second argument is simply an integer. 83
It's a whole number, and that number is basically the column of this n-dimensional array that we want 84
to look at. 85
All right. 86
So we've got a data frame, X_including_const, and say we want to calculate the variance 87
inflation factor for the first column in this data frame. 88
How would we do it? 89
Well I'll come down here 90
and as an argument pass in the first one which was "exog = X_incl 91
_const", comma, and then "exog_ids = 1". 92
Why did I pick index 1 and not 0? 93
Well at index 0 we've got our intercept, our constant. 94
I'm interested in our first feature which is Crime. 95
Okay so I'm going hit Shift+Enter on this and you're going to see that it doesn't work. 96
Take a look. 97
We get a type error. Scrolling down we see that the error description is "Unhashable type: 'slice'". 98
Why did we get this? The error message is actually fairly cryptic - unhashable type, but we do discover 99
that it actually has to do something with data types again. Now data types are actually usually quite 100
hidden in Python, 101
but occasionally, if you ignore them too much they come and bite you when you're not careful. 102
Remember how when we hit Shift+Tab here and we looked at what kind of argument this exog argument 103
should be? 104
We found out that it was an n-dimensional array. 105
Well let's look at the type of X_incl_const. 106
So writing "type(X_incl_const)" and hitting Shift+Enter will tell 107
us that this is not an n dimensional array. 108
This is in fact a dataframe, a pandas DataFrame. 109
So how do we get a n-dimensional array, an ndarray, from a dataframe? Pandas fortunately has a solution 110
for us. 111
All we have to do is call the values attribute, so all we have to do is put a dot after it and use the 112
values attribute. This retrieves an n-dimensional array from our dataframe. And just to prove that it 113
works, 114
I'm going to hit Shift+Enter and print out the value - 1.7. 1.7 is the variance inflation 115
factor for our crime feature. 116
Nice, but how do we calculate the variance inflation factors for all the other features? To do that, 117
we're gonna use a loop. We will need to loop across all the columns in the dataframe. 118
But suppose you don't know how many columns there are in your dataframe. 119
So as a challenge can you output the number of columns in X_incl_const? Bonus 120
points if you can do it in two ways. 121
I'll give you a few seconds to pause the video. Ready? 122
Here's the solution. 123
The first way you can find out how many columns there are in the dataframe are using the length function 124
and you would use the length function on the index of columns. 125
So as an argument you would pass in the dataframe and then put a dot after it and write columns, and 126
this will give you the number 14. 127
The second way you can do this is by using the shape attribute. 128
So writing the name of the dataframe and using ".shape" will give you the number of rows and the number 129
of columns. 130
So there's 404 rows and 14 columns. To access the number of columns directly use the 131
square brackets and then a one to return the second element from the shape attributes tuple. 132
There we go. 133
So both of these lines of Python code work as a way of retrieving the number of columns in a data frame. 134
Now, 135
we said we have to write a loop, 136
right? 137
So, as a second challenge, can you write a for loop that prints out all the variance inflation factors 138
on all the features? 139
I'll give you a few seconds to pause the video and try this out. 140
Ready? 141
Here's the solution. 142
You write "for i in range()", and the range would be the number of columns in the dataframe, so I could have 143
"X_incl_const.shape[1]" and then colon, and then in the body of 144
my for loop, I have a print statement with a function call to variance_inflation_factor, so "variance_inflation 145
_factor()", 146
and then my first argument is going to be, as before, "X_incl_const.values, 147
exog_ 148
idx = ", 149
And then what? Now let's apply the iterator of the loop and that's i. 150
There we go. 151
That's it. 152
That's the solution. 153
So when the loop finishes we can print "All done!" and let's run it see what we get. Tada! 154
These are all the variance inflation factors calculated for every single feature in our dataframe. 155
Now printing all this stuff out is very well and good, but what if we wanted to store it in a list? 156
What if we wanted to store all these things in a variable? 157
So I'm going to copy the code I have above, 158
come down here, add a few more cells, paste it in, delete 159
my challenge comment, 160
delete this print statement and then I'm going to create an empty list here, "vif = []", 161
gonna put a comment here, 162
this is an empty list, and then inside my list I'm going to use the append method, 163
so "vif.append()" and what will it append? It will append the variance inflation factor that is calculated 164
as part of the loop. And when we're all done, we can print our variance inflation factor list. 165
Let's see what we get. 166
So that worked really well. 167
We get all the same results printed out as before except now we're storing them in a variable. 168
But I tell you what, I'm going to show you a different kind of Python syntax to accomplish the same thing. 169
So I'm going to copy, this paste it in here. 170
What I'm going to do now is I'm going to run this loop inside these square brackets, and the way I would do 171
this is I would move this part here which is the body of the loop, 172
if you think about it, and then afterwards I will add the code for the loop, namely this part. So I can 173
copy that, paste it in here and note the square bracket at the end, 174
then I'm going to delete this part here and I will split this one line over two lines, just hitting Enter 175
here. 176
That way you can see a lot better what's going on. 177
Mm hmm. 178
Let me run this just to prove that it works. 179
Here you go. 180
So this is interesting, right? 181
We're used to seeing loops like this, but you could also use this syntax here to run the loop inside 182
the square brackets to populate your list. 183
Now, I don't really like this formatting here, so I'm going to write some code to add a dataframe at 184
the end using the dictionary notation again. 185
And this dictionary is gonna have "'coef_name':" and then the names are going to be 186
from my training dataset "X_incl_const.columns". 187
So this is gonna be my list of names, and I'm going to put comma, hit Enter to go down to a new line. For my 188
second column in my dataframe, I'll have "vif" as the column name and then I'll put the variants inflation 189
factor list that we've calculated using our loop right afterwards. 190
Now let's see what we get. We get something like this. 191
We get the feature names here and we get the variance inflation factor for each feature next to it. 192
Now there's a lot of numbers after the decimal point with the variance inflation factor. 193
So let's round it and then we use a rounding function from numpy, "np.round()", I'll have 194
the variance inflation factor inside for the first argument and for the second argument I'm going to put 195
2 as how many numbers I want after the decimal point. 196
So I want two numbers after the decimal point. I'll put a closing parentheses here and hit Shift+Enter to refresh 197
my cell. 198
And this is what we get. Now the variance inflation factors are a lot easier to read. 199
So how do we interpret this output? 200
Well, similar to how we did things with the p-value, we compare these numbers to a threshold and for the 201
variance inflation factors, the scientific consensus seems to be that the threshold is around 10, meaning 202
any feature that has a VIF over 10 would be considered problematic and would need closer inspection. 203
But looking at our list here, we can see that all the numbers are below 10 and this suggests that we 204
don't have to worry about multicollinearity. 205
Now some academics are a little bit more conservative regarding their threshold and they think that 206
a cutoff of 5 is better. 207
But given the fact that our results for the coefficients make a lot of sense and we are below the threshold 208
of 10 for the variance inflation factor, I'm really not too worried about having a multicollinearity 209
problem. Now, 210
as I've discussed earlier this housing dataset actually came from a research paper where two academics 211
were looking for the demand for clean air in Boston. 212
They were trying to put a value to how much people are willing to pay to live somewhere with good air 213
quality. 214
And I was curious if these two academics tested for multicollinearity as well 215
and if they mentioned it in their paper, and what I found was that they indeed discussed this problem 216
in the article's footnotes. 217
The problem that the researchers faced was a little different from us, because the researchers were actually 218
trying to estimate the value of having clean air and what they had was more than one measure of pollution. 219
For us in our data set we have one variable, NOX, which measures nitric oxide. But originally this 220
wasn't the only pollution factor that was tested in this research paper. 221
They actually had two different pollution measures and what they found was they actually had a 222
multicollinearity problem when they included both of them. 223
In other words, one of the pollution features was redundant, and they ended up removing it from their 224
regression model. 225
So I think this is actually quite an interesting solution to the problem 226
if you're encountering it in your own research - removing an unnecessary feature is a perfectly valid 227
way to modify the model and it's something we're going to look at in the next lesson, because in the 228
next lesson we're going to be making small tweaks to our regression model and try and simplify it and 229
at the same time we get a chance to check for that final symptom of multicollinearity, namely seeing 230
if our coefficient estimates change dramatically when we tweak the model. 231
Now it's getting late over here aand I really need to get another coffee, but I've had so much coffee 232
already so I think we'll have to maybe dilute my coffee with decaf. 233
That should work, right? 234
Anyhow, I'll see in the next lesson.