Learn how to compare and choose between linear regression models.

Choosing a Linear Regression Model

So far, we've used R-squared, adjusted R-squared, and an F-test to compare models. These criteria are most useful for finding a model that best fits an observed set of data. They are often used when our goal is interpreting a model to understand relationships between variables.

If our goal is to choose the best model for making predictions for new/unobserved data, we may want to use a likelihood based criteria instead.

*Log-likelihood* of a linear regression model essentially measures the probability of observing our data given a particular model. Higher log-likelihood is better. 

For example, we can compare two models based on log likelihood as follows:

```python
model1 = sm.OLS.from_formula('rent ~ bedrooms + size_sqft + min_to_subway', data=rentals).fit()

model2 = sm.OLS.from_formula('rent ~ bathrooms + building_age_yrs + borough', data=rentals).fit()

print(model1.llf) #Output: -44282.327
print(model2.llf) #Output: -44740.623
```

Because model 1 has a higher log-likelihood (a smaller negative number is larger), we would choose model 1 over model 2.

Log-Likelihood

In this lesson, we'll discuss some of the ways we can compare and choose linear regression models using a variety of different methods.

For example, suppose that we work at a bike rental company and have a dataset where each row represents a unique day of business (including info about the number of bikes that were rented, the weather, the day of the week, etc.). We might want to use this dataset to predict how many bikes will be rented tomorrow (so we can plan). Alternatively, we might want to use the dataset to understand which factors are most predictive of bike usage.

For either goal, we can use a linear regression model. The problem is: there are many different models that we could create. How do we know which one to use?

This lesson will focus on some common ways to compare and choose a linear model, both for prediction and data analysis.

Introduction

*R-squared* is one of the most common metrics to evaluate linear regression models. We can interpret R-squared as the proportion of variation in an outcome variable that is explained by a linear regression model. More explained variation is generally better.

For example, suppose we have a dataset containing information about apartment rentals for NYC apartments. We can build two different models to predict rental price and print out the R-Squared for each model as follows:

```python
# Create and fit the first model to predict rent
model1 = sm.OLS.from_formula('rent ~ bedrooms + size_sqft + min_to_subway', data=rentals).fit()

# Create and fit the second model
model2 = sm.OLS.from_formula('rent ~ bathrooms + building_age_yrs + borough', data=rentals).fit()

# Print out R-squared for both models
print(model1.rsquared) #Output: 0.664
print(model2.rsquared) #Output: 0.596
```
This tells us that the first model (using bedrooms, square-footage, and minutes to the subway) explains about 66.4% of the variation in rental prices, whereas the second model only explains about 59.6% of the variation. This would lead us to choose the first model over the second.

R-Squared

While R-squared is useful for comparing models with different sets of predictors, we saw that it could lead to overfitting when choosing between nested models. 

To address this issue, we can instead use adjusted R-squared, which gives a small penalty for each additional predictor in a model. For example:

```python
model1 = sm.OLS.from_formula('rent ~ bedrooms + size_sqft + borough', data=rentals).fit()

model2 = sm.OLS.from_formula('rent ~ bedrooms + size_sqft + borough + has_doorman', data=rentals).fit()

print(model1.rsquared) #Output: 0.72761
print(model2.rsquared) #Output: 0.72765

print(model1.rsquared_adj) #Output: 0.72739
print(model2.rsquared_adj) #Output: 0.72738
```
  
Note that the second model (with an additional predictor) has a slightly larger R-squared, but a slightly smaller adjusted R-squared, compared to the first model. Based on the adjusted R-squared, we would choose the smaller model.

Adjusted R-Squared

Another way of choosing a model to make predictions for new data (also called out-of-sample prediction) is by using *training* and *test* datasets. The idea is that we only use PART of our data to fit the model, then we see how well the model performs in predicting the outcome of interest for the rest of our data. The process is as follows: 

* First, we split our data into two subsets: a training set and a test set. Often, the training set is a larger proportion of the data.
* In other Python libraries, there are built-in functions to split a dataframe, but for the sake of understanding, we'll do it explicitly here by randomly sampling from a list of row indices:

```python
import numpy as np

# Create a list of indices
indices = range(len(rentals))

# Determine the size of the training set (s)
s = int(0.8*len(indices))

# Randomly select 80% of the indices
train_ind = np.random.choice(indices, size = s, replace = False)

# Create a list of the remaining 20% of indices
test_ind = list(set(indices) - set(train_ind))

# Split the data into the training and test sets
rentals_train = rentals.iloc[train_ind]
rentals_test = rentals.iloc[test_ind]
```

* Next, we fit the models we want to compare using the training set data only: 

```python
model1 = sm.OLS.from_formula('rent ~ bedrooms + bathrooms + size_sqft + min_to_subway + floor', data=rentals_train).fit()

model2 = sm.OLS.from_formula('rent ~ bedrooms + bathrooms + size_sqft + min_to_subway + floor + borough', data=rentals_train).fit()
```

* Then, we use those models to predict the rental price for the apartments in the test set:

```python
fitted1 = model1.predict(rentals_test)
fitted2 = model2.predict(rentals_test)
```

* Finally, we can compare the predicted rents to the true rents in the test set and use a metric to determine how well each model performed.
* In this example, we'll use a metric called *predictive root mean squared error (PRMSE)*, which is exactly what the name sounds like: the square root of the mean squared difference between predicted and true values of the outcome variable. A smaller PRMSE means that the model performed better (the predicted values were more similar to the true values):

```python
true = rentals_test.rent
prmse1 = np.mean((true-fitted1)**2)**.5
prmse2 = np.mean((true-fitted2)**2)**.5
print(prmse1) #output: 1326.258
print(prmse2) #output: 1224.269
```

Based on this metric, we would choose the second model over the first one because it has a smaller PRMSE.

Training and Test Sets

Similarly to R-squared, log-likelihood only increases as we add more predictors to a model. In the same way that adjusted R-squared penalizes R-squared for more predictors, there are criteria that penalize the log-likelihood for more predictors.

The two most commonly used are *Akaike information criterion (AIC)* and *Bayesian information criterion (BIC)*. Both AIC and BIC use negative log-likelihood, so we actually want the model with the LOWEST AIC or BIC.

AIC and BIC are similar, but BIC gives a bigger penalty for each additional predictor, so it is used for finding the best "simple" model. This is useful because it makes the model more interpretable. For example:

```python
model1 = sm.OLS.from_formula('rent ~ bedrooms + size_sqft + borough', data=rentals).fit()

model2 = sm.OLS.from_formula('rent ~ bedrooms + size_sqft + borough + has_doorman', data=rentals).fit()

print(model1.llf) #Output: -43756.418
print(model2.llf) #Output: -43756.017

print(model1.aic) #Output: 87522.837
print(model2.aic) #Output: 87524.034

print(model1.bic) #Output: 87555.423
print(model2.bic) #Output: 87563.137
```

We see that the log-likelihood for model 2 is slightly larger (better), but the AIC for model 2 is slightly larger (worse), and BIC even more so. Both AIC and BIC would lead us to choose model 1, whereas log-likelihood would lead us to choose model 2.

AIC and BIC

In the previous exercises, we compared nested models based on adjusted R-squared.

Another way to compare nested models is by using a hypothesis test called an F-test. Suppose we want to compare the following two models:

```python
model1 = sm.OLS.from_formula('rent ~ bedrooms + size_sqft', data=rentals).fit()

model2 = sm.OLS.from_formula('rent ~ bedrooms + size_sqft + building_age_yrs + has_elevator', data=rentals).fit()
```
Note that the second model has two more predictors than the first (`building_age_yrs` and `has_elevator`). For an F-test comparing these two models:

- The *null hypothesis* is that the coefficients on `building_age_yrs` and `has_elevator` are equal to zero (they are not useful in explaining the observed variation in rent).
- The *alternative hypothesis* is that least one of the coefficients is non-zero.

We can run the test in Python as follows:

```python
from statsmodels.stats.anova import anova_lm
anova_results = anova_lm(model1, model2)
print(anova_results)
```

Output: 

|   | df_resid | ssr     | df_diff | ss_diff | F     | Pr(>F)  |
|---|----------|---------|---------|---------|-------|---------|
| 0 | 4997.0   | 1.4e+10 | 0.0     | NaN     | NaN   | NaN     |
| 1 | 4995.0   | 1.4e+10 | 2.0     | 9.2e+08 | 170.9 | 1.6e-72 |

<br>

The p-value (`1.6e-72`, which is equal to .00000..[72 total zeros]..16) is located in the bottom right corner of this table. The column name `Pr(>F)` means "the probability of observing an F statistic greater than observed (170.9) if the null hypothesis is true".

Using a significance threshold of 0.05, the p-value is below the threshold. Therefore, we would conclude that either (or both) of the coefficients on `building_age_yrs` and `has_elevator` is non-zero. Thus, including at least one of these two predictors significantly improves the model.

This would lead us to choose `model2` over `model1`. After running this test, we might also want to separately compare a model with `building_age_yrs` and a model with `has_elevator` to see whether both are necessary. We could do this with separate F-tests or adjusted R-squared.

F-Tests

Let's again suppose that we want to use the StreetEasy data to predict rental prices in NYC. We have the following two models that we want to compare:

```python
# Fit model 1
model1 = sm.OLS.from_formula('rent ~ bedrooms + size_sqft + min_to_subway', data=rentals).fit()

# Fit model 2
model2 = sm.OLS.from_formula('rent ~ bedrooms + size_sqft + min_to_subway + borough', data=rentals).fit()

# Print out R-squared for both models
print(model1.rsquared) # Output: 0.664
print(model2.rsquared) # Output: 0.728
```
Note that these models both use `bedrooms`, `size_sqft`, and `min_to_subway` as predictors; but `model2` uses `borough` as well. Because all of the predictors in `model1` are also contained in `model2`, these are called *nested models*.

It turns out that larger nested models will ALWAYS have higher R-squared than their smaller counterparts. However, adding a lot of additional predictors can lead to a different issue: over-fitting. To understand the intuition behind why overfitting is problematic, consider the following plot of rental prices vs. number of bathrooms. We can perfectly predict each datapoint if we fit the zig-zagging line shown below:

![plot showing rent vs. bedrooms for 9 apartments. The points are all connected by a zig-zagging dotted line.](https://static-assets.codecademy.com/skillpaths/master-stats-ii/choosing-linear-regression-model/overfitting.svg)

However, imagine that we collect a new sample of apartments in NYC and record the number of bathrooms in each. Then, suppose we want to use our model to predict rental prices. Even if the overall relationship between bathrooms and rent is the same in our new data, the exact values will be slightly different. Predictions based on the zig-zag line may be less accurate because the model was so heavily influenced by the quirks of the data we originally collected. A straight line through the middle of the points is actually more useful.

R-Squared and the Dangers of Overfitting

Congratulations! In this lesson, you've learned a number of different methods for model comparison:

* For choosing a model that best represents the data we have:
  * R-squared
  * Adjusted R-squared
  * F-test
* For choosing a model for accurate out-of-sample prediction:
  * Log likelihood
  * AIC/BIC
  * Training/test sets

Note that we've covered many different methods for choosing a model and they don't always agree. In order to choose a method, it's important to consider your ultimate goal (analysis vs. prediction) and what you want to prioritize (simplicity and interpretability vs. accuracy)