<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Introduction to Linear Regression
Adapted from Chapter 3 of An Introduction to Statistical Learning
Predictive modeling, using a data samples to make predictions about unobserved or future events, is a common data analytics task. Predictive modeling is considered to be a form of machine learning.
Linear regression is a technique for predicting a response/dependent variable based on one or more explanatory/independent variables, or features. The term "linear" refers to the fact that the method models data as a linear combination of explanatory variables.  
Linear regression, in its simplest form, fits a straight line to the response variable data so that the line minimizes the squared differences (also called errors or residuals) between the actual obbserved response and the predicted point on the line. Since linear regression fits the observed data with a line, it is most effective when the response and the explanatory variable do have a linear relationship.
Motivating Example
Step1:  The features?
- TV
Step2:  There are 200 observations, corresponding to 200 markets.
We can try to discover if there is any relationship between the money spend on a specific type of ad, in a given market, and the sales in that market by plotting the sales figures against each category of advertising expenditure.
Step3:  Questions
How can the company selling the product decide on how to spend its advertising money in the future? We first need to answer the following question
Step4:  Interpreting Model Coefficients
Q
Step5:  The predicted Sales in that market are of 9.409444 * 1000 =~ 9409 widgets
Using Statsmodels
Step6:  Plotting the Least Squares Line
Let's make predictions for the smallest and largest observed values of money spent on TV ads, and then use the predicted values to plot the least squares line
Step7:  Confidence in Linear Regression Models
Q
Step8:  Since we only have a single sample of data, and not the entire population the "true" value of the regression coefficient is either within this interval or it isn't, but there is no way to actually know. 
We estimate the regression coefficient using the data we have, and then we characterize the uncertainty about that estimate by giving a confidence interval, an interval that will "probably" contain the  value coefficient. Note that there is no probability associated with the true value of the regression coefficient being in the given confidence interval!
Also note that using 95% confidence intervals is simply a convention. One can create 90% confidence intervals (narrower intervals), 99% confidence intervals (wider intervals), etc.
Hypothesis Testing and p-values
Closely related to confidence intervals is hypothesis testing. Generally speaking, you start with a null hypothesis and an alternative hypothesis (that is opposite the null). Then, you check whether the data supports rejecting the null hypothesis or failing to reject the null hypothesis.
(Note that "failing to reject" the null is not the same as "accepting" the null hypothesis. The alternative hypothesis may indeed be true, except that you just don't have enough data to show that.)
As it relates to model coefficients, here is the conventional hypothesis test
Step9:  If the 95% confidence interval includes zero, the p-value for that coefficient will be greater than 0.05. If the 95% confidence interval does not include zero, the p-value will be less than 0.05. Thus, a p-value less than 0.05 is one way to decide whether there is likely a relationship between the feature and the response. (Again, using 0.05 as the cutoff is just a convention.)
In this case, the p-value for TV is far less than 0.05, and so we believe that there is a relationship between TV ads and Sales.
Note that we generally ignore the p-value for the intercept.
How Well Does the Model Fit the data?
The most common way to evaluate the overall fit of a linear model to the available data is by calculating the R-squared (a.k.a, "coefficient of determination") value. 
R-squared has several interpretations
Step10:  Is that a "good" R-squared value? One cannot generally assess that. What a "good" R-squared value is depends on the domain and therefore R-squared is most useful as a tool for comparing different models.
Multiple Linear Regression
Simple linear regression can be extended to include multiple explanatory variables
Step11:  How do we interpret the coefficients? For a given amount of Radio and Newspaper ad spending, an increase of a unit ($1000 dollars) in TV ad spending is associated with an increase in Sales of 45.765 widgets.
Other information is available in the model summary output
Step12:  TV and Radio have significant p-values, whereas Newspaper does not. Thus we reject the null hypothesis for TV and Radio (that there is no association between those features and Sales), and fail to reject the null hypothesis for Newspaper.
TV and Radio ad spending are both positively associated with Sales, whereas Newspaper ad spending is slightly negatively associated with Sales. 
This model has a higher R-squared (0.897) than the previous model, which means that this model provides a better fit to the data than a model that only includes TV.
Feature Selection
How do I decide which features to include in a linear model? 
- Try different models and check whether the R-squared value goes up when you add new predictors.
What are the drawbacks to this approach?
- Linear models rely upon a lot of assumptions (such as the predictors/features being independent), and if those assumptions are violated (which they usually are), R-squared are less reliable.
- R-squared is susceptible to overfitting, and thus there is no guarantee that a model with a high R-squared value will generalize well to new data. For example
Step13:  R-squared will always increase as you add more features to the model, even if they are unrelated to the response. Thus, selecting the model with the highest R-squared is not a reliable approach for choosing the best linear model.
There is alternative to R-squared called adjusted R-squared that penalizes model complexity (to control for overfitting), but this approach has its own set of issues.
Is there a better approach to feature selection? Cross-validation, which provides a more reliable estimate of out-of-sample error, and thus is better at choosing which model will better generalize to out-of-sample data. Cross-validation can be applied to any type of model, not just linear models.
Linear Regression in scikit-learn
The work done using Statsmodels can also be using scikit-learn
Step14:  Handling Categorical Predictors with Two Categories
What if one of the predictors was categorical?
Let's create a new feature called Size, and randomly assign observations to be small or large
Step15:  When using scikit-learn, we need to represent all data numerically. For example, if the feature we want to represent has only two categories, we create a dummy variable that represents the categories as a binary value
Step16:  The multiple linear regression including the IsLarge predictor
Step17:  How do we interpret the coefficient of IsLarge? For a given amount of TV/Radio/Newspaper ad spending, a large market is associated with an average increase in Sales of 51.55 widgets (as compared to sales in a Small market).
If we reverse the 0/1 encoding and created the feature 'IsSmall', the coefficient would be the same in absolute value, but negative instead of positive. All that changes is the interpretation of the coefficient.
Handling Categorical Predictors with More than Two Categories
Let's create a new feature called Area, and randomly assign observations to be rural, suburban, or urban
Step18:  We have to represent Area numerically, but an encoding such as 0=rural, 1=suburban, 2=urban would not work because that would imply that there is an ordered relationship between suburban and urban. Instead, we can create another dummy variable.
Step19:  rural is coded as Area_suburban=0 and Area_urban=0
suburban is coded as Area_suburban=1 and Area_urban=0
urban is coded as Area_suburban=0 and Area_urban=1
Only two dummies are needed to captures all of the information about the Area feature.(In general, for a categorical feature with k levels, we create k-1 dummy variables.)
Let's include the two new dummy variables in the model
<ASSISTANT_TASK:>
Python Code:
# imports
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline
# read data into a DataFrame
data = pd.read_csv('http://www-bcf.usc.edu/~gareth/ISL/Advertising.csv', index_col=0)
data.head()
Explanation: Introduction to Linear Regression
Adapted from Chapter 3 of An Introduction to Statistical Learning
Predictive modeling, using a data samples to make predictions about unobserved or future events, is a common data analytics task. Predictive modeling is considered to be a form of machine learning.
Linear regression is a technique for predicting a response/dependent variable based on one or more explanatory/independent variables, or features. The term "linear" refers to the fact that the method models data as a linear combination of explanatory variables.  
Linear regression, in its simplest form, fits a straight line to the response variable data so that the line minimizes the squared differences (also called errors or residuals) between the actual obbserved response and the predicted point on the line. Since linear regression fits the observed data with a line, it is most effective when the response and the explanatory variable do have a linear relationship.
Motivating Example: Advertising Data
Let us look at data depicting the money(in thousands of dollars) spent on TV, Radio and newspaper ads for a product in a given market, as well as the corresponding sales figures.
End of explanation
# print the size of the DataFrame object, i.e., the size of the dataset
data.shape
Explanation: The features?
- TV: advertising dollars spent on TV (for a single product, in a given market) 
- Radio: advertising dollars spent on Radio (for a single product, in a given market) 
- Newspaper: advertising dollars spent on Newspaper (for a single product, in a given market) 
What is the response?
- Sales: sales of a single product in a given market (in thousands of widgets)
End of explanation
fig, axs = plt.subplots(1, 3, sharey=True)
data.plot(kind='scatter', x='TV', y='Sales', ax=axs[0], figsize=(16, 8))
data.plot(kind='scatter', x='Radio', y='Sales', ax=axs[1])
data.plot(kind='scatter', x='Newspaper', y='Sales', ax=axs[2])
Explanation: There are 200 observations, corresponding to 200 markets.
We can try to discover if there is any relationship between the money spend on a specific type of ad, in a given market, and the sales in that market by plotting the sales figures against each category of advertising expenditure.
End of explanation
import statsmodels.formula.api as sf
#create a model with Sales as dependent variable and TV as explanatory variable
model = sf.ols('Sales ~ TV', data)
#fit the model to the data 
fitted_model = model.fit()
# print the coefficients
print(fitted_model.params)
Explanation: Questions
How can the company selling the product decide on how to spend its advertising money in the future? We first need to answer the following question: "Based on this data, does there apear to be a relationship between ads and sales?"
If yes, 
1. Which ad types contribute to sales?
2. How strong is the relationship between each ad type and sales?
4. What is the effect of each ad type of sales?
5. Given ad spending in a particular market, can sales be predicted?
We will use Linear Regression to try and asnwer these questions.
Simple Linear Regression
Simple linear regression is an approach for modeling the relatrionship between a dependent variable (a "response") and an explanatory variable, also known as a "predictor" or "feature". The relationship is modeled as a linear function $y = \beta_0 + \beta_1x$ whose parameters are estimated from the available data.
In the equation above:
- $y$ is called the response, regressand, endogenous variable, dependent variable, etc.
- $x$ is the feature, regressor, exogenous variable, explanatory variables, predictor, etc.
- $\beta_0$ is known as the intercept
- $\beta_1$ is the regression coefficient, effect, etc. 
Together, $\beta_0$ and $\beta_1$ are called paramaters, model/regression coefficients, or effects. To create a model, we must discover/learn/estimate the values of these coefficients. 
Estimating/Learning Model/Regression Coefficients
Regression coefficients are estimated using a variety of methods. The least squares method, which finds the line which minimizes the sum of squared residuals (or "sum of squared errors") is among the most oftenly used.
In the pictures below:
- The blue dots are the observed values of x and y.
- The red line is the least squares line.
- The residuals are the distances between the observed values and the least squares line.
$\beta_0$ is the intercept of the least squares line (the value of $y$ when $x$=0)
$\beta_1$ is the slope of the least squares line, i.e. the ratio of the vertical change (in $y$) and the horizontal change (in $x$).
We can use the  statsmodels package to estimate the model coefficients for the advertising data:
End of explanation
7.032594 + 0.047537*50
Explanation: Interpreting Model Coefficients
Q: How do we interpret the coefficient ($\beta_1$) of the explanatory variable "TV"?
A: A unit (a thousand dollars) increase in TV ad spending is associated with a 0.047537 unit (a thousand widgets) increase in Sales, i.e., an additional $1000 spent on TV ads is associated with an increase in sales of ~47.5 widgets.
Note that it is, in general, possible to have a negative effect, e.g., an increase in TV ad spending to be associated with a decrease in sales. $\beta_1$ would be negative in this case.
Using the Model for Prediction
Can we use the model we develop to guide advertising spending decisions? For example, if the company spends $50,000 on TV advertising in a new market, what would the model predict for the sales in that market?
$$y = \beta_0 + \beta_1x$$
$$y = 7.032594 + 0.047537 \times 50$$
End of explanation
# create a DataFrame to use with the Statsmodels formula interface
New_TV_spending = pd.DataFrame({'TV': [50]})
#check the newly created DataFrame
New_TV_spending.head()
# use the model created above to predict the sales to be generated by the new TV ad money
sales = fitted_model.predict(New_TV_spending)
print(sales)
Explanation: The predicted Sales in that market are of 9.409444 * 1000 =~ 9409 widgets
Using Statsmodels:
End of explanation
# create a DataFrame with the minimum and maximum values of TV ad money
New_TV_money = pd.DataFrame({'TV': [data.TV.min(), data.TV.max()]})
print(New_TV_money.head())
# make predictions for those x values and store them
sales_predictions = fitted_model.predict(New_TV_money)
print(sales_predictions)
# plot the observed data
data.plot(kind='scatter', x='TV', y='Sales')
# plot the least squares line
plt.plot(New_TV_money, sales_predictions, c='red', linewidth=2)
Explanation: Plotting the Least Squares Line
Let's make predictions for the smallest and largest observed values of money spent on TV ads, and then use the predicted values to plot the least squares line:
End of explanation
# print the confidence intervals for the model coefficients
print(fitted_model.conf_int())
Explanation: Confidence in Linear Regression Models
Q: Is linear regression a high bias/low variance model, or a low variance/high bias model?
A: High bias/low variance. Under repeated sampling, the line will stay roughly in the same place (low variance), but the average of those models won't do a great job capturing the true relationship (high bias). (A low variance is a useful characteristic when limited training data is available.)
We can use Statsmodels to calculate 95% confidence intervals for the model coefficients, which are interpreted as follows: If the population from which this sample was drawn was sampled 100 times, approximately 95 of those confidence intervals would contain the "true" coefficient.
End of explanation
# print the p-values for the model coefficients
fitted_model.pvalues
Explanation: Since we only have a single sample of data, and not the entire population the "true" value of the regression coefficient is either within this interval or it isn't, but there is no way to actually know. 
We estimate the regression coefficient using the data we have, and then we characterize the uncertainty about that estimate by giving a confidence interval, an interval that will "probably" contain the  value coefficient. Note that there is no probability associated with the true value of the regression coefficient being in the given confidence interval!
Also note that using 95% confidence intervals is simply a convention. One can create 90% confidence intervals (narrower intervals), 99% confidence intervals (wider intervals), etc.
Hypothesis Testing and p-values
Closely related to confidence intervals is hypothesis testing. Generally speaking, you start with a null hypothesis and an alternative hypothesis (that is opposite the null). Then, you check whether the data supports rejecting the null hypothesis or failing to reject the null hypothesis.
(Note that "failing to reject" the null is not the same as "accepting" the null hypothesis. The alternative hypothesis may indeed be true, except that you just don't have enough data to show that.)
As it relates to model coefficients, here is the conventional hypothesis test:
- null hypothesis: There is no relationship between TV ads and Sales (and thus $\beta_1$ equals zero)
- alternative hypothesis: There is a relationship between TV ads and Sales (and thus $\beta_1$ is not equal to zero)
How do we test this hypothesis? Intuitively, we reject the null (and thus believe the alternative) if the 95% confidence interval does not include zero. Conversely, the p-value represents the probability that the coefficient is actually zero:
End of explanation
# print the R-squared value for the model
fitted_model.rsquared
Explanation: If the 95% confidence interval includes zero, the p-value for that coefficient will be greater than 0.05. If the 95% confidence interval does not include zero, the p-value will be less than 0.05. Thus, a p-value less than 0.05 is one way to decide whether there is likely a relationship between the feature and the response. (Again, using 0.05 as the cutoff is just a convention.)
In this case, the p-value for TV is far less than 0.05, and so we believe that there is a relationship between TV ads and Sales.
Note that we generally ignore the p-value for the intercept.
How Well Does the Model Fit the data?
The most common way to evaluate the overall fit of a linear model to the available data is by calculating the R-squared (a.k.a, "coefficient of determination") value. 
R-squared has several interpretations:
(1) R-squared ×100 percent of the variation in the dependent variable ($y$) is reduced by taking into account predictor $x$
(2) R-squared is the proportion of variance in the observed data that is "explained" by the model.
R-squared is between 0 and 1, and, generally speaking, higher is considered to be better because more variance is accounted for ("explained") by the model.
Note, however, that R-squared does not indicate whether a regression model is actually good. You can have a low R-squared value for a good model, or a high R-squared value for a model that does not fit the data!
One should evaluate the adequacy of a model by looking at R-squared values as well as residual (i.e., observed value - fitted value) plots, other model statistics, and subject area knowledge.
The R-squared value for our simple linear regression model is:
End of explanation
# create a model with all three features
multi_model = sf.ols(formula='Sales ~ TV + Radio + Newspaper', data=data)
fitted_multi_model = multi_model.fit()
# print the coefficients
print(fitted_multi_model.params)
Explanation: Is that a "good" R-squared value? One cannot generally assess that. What a "good" R-squared value is depends on the domain and therefore R-squared is most useful as a tool for comparing different models.
Multiple Linear Regression
Simple linear regression can be extended to include multiple explanatory variables:
$y = \beta_0 + \beta_1x_1 + ... + \beta_nx_n$
Each $x$ represents a different predictor/feature, and each predictor has its own coefficient. In our case:
$y = \beta_0 + \beta_1 \times TV + \beta_2 \times Radio + \beta_3 \times Newspaper$
Let's use Statsmodels to estimate these coefficients:
End of explanation
# print a summary of the fitted model
fitted_multi_model.summary()
Explanation: How do we interpret the coefficients? For a given amount of Radio and Newspaper ad spending, an increase of a unit ($1000 dollars) in TV ad spending is associated with an increase in Sales of 45.765 widgets.
Other information is available in the model summary output:
End of explanation
# only include TV and Radio in the model
model1 = sf.ols(formula='Sales ~ TV + Radio', data=data).fit()
print(model1.rsquared)
# add Newspaper to the model (which we believe has no association with Sales)
model2 = sf.ols(formula='Sales ~ TV + Radio + Newspaper', data=data).fit()
print(model2.rsquared)
Explanation: TV and Radio have significant p-values, whereas Newspaper does not. Thus we reject the null hypothesis for TV and Radio (that there is no association between those features and Sales), and fail to reject the null hypothesis for Newspaper.
TV and Radio ad spending are both positively associated with Sales, whereas Newspaper ad spending is slightly negatively associated with Sales. 
This model has a higher R-squared (0.897) than the previous model, which means that this model provides a better fit to the data than a model that only includes TV.
Feature Selection
How do I decide which features to include in a linear model? 
- Try different models and check whether the R-squared value goes up when you add new predictors.
What are the drawbacks to this approach?
- Linear models rely upon a lot of assumptions (such as the predictors/features being independent), and if those assumptions are violated (which they usually are), R-squared are less reliable.
- R-squared is susceptible to overfitting, and thus there is no guarantee that a model with a high R-squared value will generalize well to new data. For example:
End of explanation
# create a DataFrame
feature_cols = ['TV', 'Radio', 'Newspaper']
X = data[feature_cols]
y = data.Sales
from sklearn.linear_model import LinearRegression
lm = LinearRegression()
lm.fit(X, y)
# print intercept and coefficients
print(lm.intercept_)
print(lm.coef_)
# pair the feature names with the coefficients
print(zip(feature_cols, lm.coef_))
# predict for a new observation
lm.predict([[100, 25, 25]])
# calculate the R-squared
lm.score(X, y)
Explanation: R-squared will always increase as you add more features to the model, even if they are unrelated to the response. Thus, selecting the model with the highest R-squared is not a reliable approach for choosing the best linear model.
There is alternative to R-squared called adjusted R-squared that penalizes model complexity (to control for overfitting), but this approach has its own set of issues.
Is there a better approach to feature selection? Cross-validation, which provides a more reliable estimate of out-of-sample error, and thus is better at choosing which model will better generalize to out-of-sample data. Cross-validation can be applied to any type of model, not just linear models.
Linear Regression in scikit-learn
The work done using Statsmodels can also be using scikit-learn:
End of explanation
import numpy as np
# create a Series of booleans in which roughly half are True
#generate len(data) numbers between 0 and 1
numbers = np.random.rand(len(data))
#create and index of 0s and 1s by based on whether the corresponding random number
#is greater than 0.5. 
index_for_large = (numbers > 0.5)
#create a new data column called Size and set its values to 'small'
data['Size'] = 'small'
# change the values of Size to 'large' whenever the corresponding value of the index is 1 
data.loc[index_for_large, 'Size'] = 'large'
data.head()
Explanation: Handling Categorical Predictors with Two Categories
What if one of the predictors was categorical?
Let's create a new feature called Size, and randomly assign observations to be small or large:
End of explanation
# create a new Series called IsLarge
data['IsLarge'] = data.Size.map({'small':0, 'large':1})
data.head()
Explanation: When using scikit-learn, we need to represent all data numerically. For example, if the feature we want to represent has only two categories, we create a dummy variable that represents the categories as a binary value:
End of explanation
# create X and y
feature_cols = ['TV', 'Radio', 'Newspaper', 'IsLarge']
X = data[feature_cols]
y = data.Sales
# instantiate, fit
lm = LinearRegression()
lm.fit(X, y)
# print coefficients
list(zip(feature_cols, lm.coef_))
Explanation: The multiple linear regression including the IsLarge predictor:
End of explanation
# set a seed for reproducibility
np.random.seed(123456)
# assign roughly one third of observations to each group
nums = np.random.rand(len(data))
mask_suburban = (nums > 0.33) & (nums < 0.66)
mask_urban = nums > 0.66
data['Area'] = 'rural'
data.loc[mask_suburban, 'Area'] = 'suburban'
data.loc[mask_urban, 'Area'] = 'urban'
data.head()
Explanation: How do we interpret the coefficient of IsLarge? For a given amount of TV/Radio/Newspaper ad spending, a large market is associated with an average increase in Sales of 51.55 widgets (as compared to sales in a Small market).
If we reverse the 0/1 encoding and created the feature 'IsSmall', the coefficient would be the same in absolute value, but negative instead of positive. All that changes is the interpretation of the coefficient.
Handling Categorical Predictors with More than Two Categories
Let's create a new feature called Area, and randomly assign observations to be rural, suburban, or urban:
End of explanation
# create three dummy variables using get_dummies, then exclude the first dummy column
area_dummies = pd.get_dummies(data.Area, prefix='Area').iloc[:, 1:]
# concatenate the dummy variable columns onto the original DataFrame (axis=0 means rows, axis=1 means columns)
data = pd.concat([data, area_dummies], axis=1)
data.head()
Explanation: We have to represent Area numerically, but an encoding such as 0=rural, 1=suburban, 2=urban would not work because that would imply that there is an ordered relationship between suburban and urban. Instead, we can create another dummy variable.
End of explanation
# read data into a DataFrame
#data = pd.read_csv('http://www-bcf.usc.edu/~gareth/ISL/Advertising.csv', index_col=0)
# create X and y
feature_cols = ['TV', 'Radio', 'Newspaper', 'IsLarge', 'Area_suburban', 'Area_urban']
X = data[feature_cols]
y = data.Sales
# instantiate, fit
lm = LinearRegression()
lm.fit(X, y)
# print coefficients
list(zip(feature_cols, lm.coef_))
Explanation: rural is coded as Area_suburban=0 and Area_urban=0
suburban is coded as Area_suburban=1 and Area_urban=0
urban is coded as Area_suburban=0 and Area_urban=1
Only two dummies are needed to captures all of the information about the Area feature.(In general, for a categorical feature with k levels, we create k-1 dummy variables.)
Let's include the two new dummy variables in the model:
End of explanation
<END_TASK>