grouped_test2.get_group('4wd')['price']

4      17450.0
136     7603.0
140     9233.0
141    11259.0
144     8013.0
145    11694.0
150     7898.0
151     8778.0
Name: price, dtype: float64

# ANOVA
f_val, p_val = stats.f_oneway(grouped_test2.get_group('fwd')['price'], grouped_test2.get_group('rwd')['price'], grouped_test2.get_group('4wd')['price'])  
 
print( "ANOVA results: F=", f_val, ", P =", p_val)   

ANOVA results: F= 67.9540650078 , P = 3.39454435772e-23

f_val, p_val = stats.f_oneway(grouped_test2.get_group('fwd')['price'], grouped_test2.get_group('rwd')['price'])  
 
print( "ANOVA results: F=", f_val, ", P =", p_val )

ANOVA results: F= 130.553316096 , P = 2.23553063557e-23

f_val, p_val = stats.f_oneway(grouped_test2.get_group('4wd')['price'], grouped_test2.get_group('rwd')['price'])  
   
print( "ANOVA results: F=", f_val, ", P =", p_val)   

ANOVA results: F= 8.58068136892 , P = 0.00441149221123

f_val, p_val = stats.f_oneway(grouped_test2.get_group('4wd')['price'], grouped_test2.get_group('fwd')['price'])  
 
print("ANOVA results: F=", f_val, ", P =", p_val)   

ANOVA results: F= 0.665465750252 , P = 0.416201166978

from sklearn.linear_model import LinearRegression

lm = LinearRegression()

X = df[['predictor']]
Y = df[['target']]

lm.fit(X,Y)

yhat = lm.predict(X)

df.head()

from sklearn.linear_model import LinearRegression

X = df[['highway-mpg']]
Y = df['price']

lm = LinearRegression()
lm.fit(X,Y)

LinearRegression(copy_X=True, fit_intercept=True, n_jobs=1, normalize=False)

Yhat=lm.predict(X)
Yhat[0:5]

array([ 16236.50464347,  16236.50464347,  17058.23802179,  13771.3045085 ,
        20345.17153508])

lm_eq = 'yhat = ' + str(lm.intercept_) + ' + ' + str(lm.coef_[0]) + 'x'
print(lm_eq)

yhat = 38423.3058582 + -821.733378322x

Multiple Linear Regression

MLR is used to explain the relationship between 1 continuous target and 2 or more predictors

The model will then be defined by the function

$\hat y=b_0 + b_1x_1 + b_2x_2 + b_3x_3 + ... + b_nx_n$

The resulting equation represents a surface in $n$ dimensional space that will define $y$ given each component of $x$

We can train this model just as before as follows

X = df[['predictor1','predictor2','predictor3',...,'predictorN']]
Y = df[['target']]

lm.fit(X,Y)

And then predict as before with

yhat = lm.predict(X)

Where X is in the form of the training data

The intercept and coefficiencts of the model can once again be found with lm.intercept_ and lm.coef_ respectively

We can develop a MLR model as follows

Z = df[['horsepower', 'curb-weight', 'engine-size', 'highway-mpg']]
lm.fit(Z, df['price'])
print('intercept: ' + str(lm.intercept_), '\ncoefficients: ' + str(lm.coef_))

intercept: -15806.6246263 
coefficients: [ 53.49574423   4.70770099  81.53026382  36.05748882]

Model Evaluation using Visualization

Regression Plot

Regression plots gives us a good estimate of:

• The relationship between two variables
• The strength of correlation
• The direction of the relationship

A regression plot is a combination of a scatter plot, and a linear regression line

We create a regression plot using the sns.regplot() function in seaborn

import seaborn as sns

sns.regplot(x='predictor', y='target', data=df)
plt.ylim(0,)
plt.show()

Single Linear Regression

We can visualize SLR as follows

width = 12
height = 10
plt.figure(figsize=(width, height))
sns.regplot(x="highway-mpg", y="price", data=df)
plt.ylim(0,)
plt.show()

plt.figure(figsize=(width, height))
sns.regplot(x="peak-rpm", y="price", data=df)
plt.ylim(0,)

(0, 47422.919330307624)

Residual Plot

A Residual plot represents the error between the actual and predicted value, the results should have 0 mean if the linear assumption is applicable

If the residual plot is not equally distributed around the mean of 0 throughout, we know that the linear model is not correct for the data we have

We can create a residual plot as with sns.residplot()

sns.residplot(df['predictor'], df['target'])

width = 12
height = 10
plt.figure(figsize=(width, height))
sns.residplot(df['highway-mpg'], df['price'])
plt.show()

Distribution Plot

We use this to look at the distribution of the actual vs predicted values for our model

A distribution plot can be made with sns.distplot()

ax1 = sns.distplot(df['target'], hist=False, label='Actual Value')
sns.distplot(yhat, hist=False, label='Fitted Values', ax=ax1)
plt.show()

Multiple Linear Regression

Using A Distribution Plot is more valuable when looking at MLR model performance

Y_hat = lm.predict(Z)

plt.figure(figsize=(width, height))


ax1 = sns.distplot(df['price'], hist=False, color="r", label="Actual Value")
sns.distplot(Yhat, hist=False, color="b", label="Fitted Values" , ax=ax1)


plt.title('Actual vs Fitted Values for Price')
plt.xlabel('Price (in dollars)')
plt.ylabel('Proportion of Cars')

plt.show()
plt.close()

Polynomial Regression and Pipelines

Polynomial Regression is a form of regression used to describe curvilinear relationships in which our predictor variable is not linear

• Quadratic $\hat y=b_0 + b_1x_1 + b_2x_1^2$
• Cubic $\hat y=b_0 + b_1x_1 + b_2x_1^2 + b_2x_1^3$
• Higher Order $\hat y=b_0 + b_1x_1 + b_2x_1^2 + ... + b_2x_1^n$

Picking the correct order can make our model more accurate

We can fit a polynomial to our data using np.polyfit(), and can print out the resulting model with np.polydl()

For an $n^{th}$ order polynomial we can create and view a model as follows

model = np.polyfit(X,Y,n)
print(np.polydl(f))

Visualization

Visualizing Polynomial regression plots is a bit more work but can be done with the function below

def PlotPolly(model,independent_variable,dependent_variabble, Name):
    x_new = np.linspace(15, 55, 100)
    y_new = model(x_new)

    plt.plot(independent_variable,dependent_variabble,'.', x_new, y_new, '-')
    plt.title('Polynomial Fit with Matplotlib for Price ~ Length')
    ax = plt.gca()
    fig = plt.gcf()
    plt.xlabel(Name)
    plt.ylabel('Price of Cars')

    plt.show()
    plt.close()

Next we can do the polynomial fit for our data

x = df['highway-mpg']
y = df['price']

f = np.polyfit(x, y, 3)
p = np.poly1d(f)
print(p)

        3         2
-1.557 x + 204.8 x - 8965 x + 1.379e+05

PlotPolly(p,x,y, 'highway-mpg')

Multi Dimension Polynomial Regression

Polynomial regression can also be done in multiple dimensions, for example a second order approximation would look like the following

$\hat y=b_0 + b_1x_1 + b_2x_2 + b_3x_1x_2 + b_4x_1^2 + b_5x_2^2$

If we want to do multi dimentional polynomial fits, we will need to use PolynomialFeatures from sklearn.preprocessing to preprocess our features as follows

X = df[['predictor1','predictor2','predictor3',...,'predictorN']]
Y = df[['target']]

pr = PolynomialFeatures(degree=2)
x_poly = pr.fit_transform(X, include_bias=False)

from sklearn.preprocessing import PolynomialFeatures

pr = PolynomialFeatures(degree=2)
pr

PolynomialFeatures(degree=2, include_bias=True, interaction_only=False)

Z_pr = pr.fit_transform(Z)

Before the transformation we have 4 features

Z.shape

(201, 4)

After the transformation we have 15 features

Z_pr.shape

(201, 15)

Pre-Processing

sklearn has some preprocessing functionality such as

Normalization

We can train a scaler and normalize our data based on that as follows

from sklearn.preprocessing import StandardScaler

scale = StandardScaler()
scale.fit(X[['feature1','feature2',..., 'featureN']])

X_scaled = scale.transform(X[['feature1','feature2',..., 'featureN']])

There are other normalization functions available with which we can preprocess our data

Pipelines

There are many steps to getting a prediction, such as Normalization, Polynomial Transformation, and training a Linear Regression Model

We can use a Pipeline library to help simplify the process

First we import all the libraries we will need as well as the pipeline library

from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import StandardScaler

from sklearn.pipeline import Pipeline

Then we construct our pipeline using a list of tuples defined as ('name of estimator model', ModelConstructor()) and create our pipeline object

input = [('scale', StandardScaler()),
         ('polynomial', PolynomialFeatures(degree=n),
         ...,
         ('mode', LinearRegression()]
pipe = Pipeline(Input)

We can then create our pipeline object on the data by using the pipe.train function

pipe.train(X[['feature1','feature2',..., 'featureN']], Y)
yhat = pipe.predict(X[['feature1','feature2',..., 'featureN']])

The above method will normalize the data, then perform a polynomial transform and output a prediction

from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler

Input=[('scale',StandardScaler()),
       ('polynomial', PolynomialFeatures(include_bias=False)),
       ('model',LinearRegression())]

pipe=Pipeline(Input)
pipe

Pipeline(memory=None,
     steps=[('scale', StandardScaler(copy=True, with_mean=True, with_std=True)), ('polynomial', PolynomialFeatures(degree=2, include_bias=False, interaction_only=False)), ('model', LinearRegression(copy_X=True, fit_intercept=True, n_jobs=1, normalize=False))])

pipe.fit(Z,y)

Pipeline(memory=None,
     steps=[('scale', StandardScaler(copy=True, with_mean=True, with_std=True)), ('polynomial', PolynomialFeatures(degree=2, include_bias=False, interaction_only=False)), ('model', LinearRegression(copy_X=True, fit_intercept=True, n_jobs=1, normalize=False))])

ypipe=pipe.predict(Z)
ypipe[0:4]

array([ 13102.74784201,  13102.74784201,  18225.54572197,  10390.29636555])

In-Sample Evaluation

Two important measures that determine how well a model fits a speficic dataset are

• Mean Squared Error (MSE)
• R-Squared ($R^2$)

Mean Squared Error

We simply find the difference between the average square error of our prediction when compared to our data

To get the MSE in Python we can do the following

from sklearn.metrics import mean_squared_error

mean_squared_error(X['target'],Y_predicted_sample)

R-Squared

$R^2$ is the coefficient of determination

• Measure of how close the data is to the fitted regression line
• The percentage of variation in Y that is explained by the model
• Like comparing our model to the mean of the data in approximating the data

$R^2$ is usually between 0 and 1

$R^2=1-\frac{MSE of Regression Line}{MSE of \bar y}$

We can get the $R^24$ value with lm.score()

A negative $R^2$ can be a sign of overfitting

In-Sample Evaluation of Models

We can evaluate our models with the following

#highway_mpg_fit
lm.fit(X, Y)
# Find the R^2
lm.score(X, Y)

0.49659118843391747

Yhat = lm.predict(X)
Yhat[0:4]

array([ 16236.50464347,  16236.50464347,  17058.23802179,  13771.3045085 ])

from sklearn.metrics import mean_squared_error

#mean_squared_error(Y_true, Y_predict)
mean_squared_error(df['price'], Yhat)

31635042.944639895

# fit the model 
lm.fit(Z, df['price'])
# Find the R^2
lm.score(Z, df['price'])

0.80935628065774567

Y_predict_multifit = lm.predict(Z)

mean_squared_error(df['price'], Y_predict_multifit)

11980366.87072649

Now we can check what the $R^2$ value for our model is

from sklearn.metrics import r2_score

r_squared = r2_score(y, p(x))
r_squared

0.67419466639065173

mean_squared_error(df['price'], p(x))

20474146.426361222

Prediction and Decision Making

We should use visualization, numerical evaluation, and model comparison in order to see if the model values makes sense

To compare our model to the data we can simply plot the output of our model over the range of our data

new_input=np.arange(1,100,1).reshape(-1,1)

lm.fit(X, Y)

LinearRegression(copy_X=True, fit_intercept=True, n_jobs=1, normalize=False)

yhat=lm.predict(new_input)
yhat[0:5]

array([ 37601.57247984,  36779.83910151,  35958.10572319,  35136.37234487,
        34314.63896655])

plt.plot(new_input,yhat)
plt.show()

Conclusion

From the results above (yes, they're a mess but it's all pretty much just from the CC Lab file) we can note the $R^2$ and MSE values are as follows

Simple Linear Regression: Using Highway-mpg as a Predictor Variable of Price.

• R-squared: 0.49659118843391759
• MSE: 3.16 x10^7

Multiple Linear Regression: Using Horsepower, Curb-weight, Engine-size, and Highway-mpg as Predictor Variables of Price.

• R-squared: 0.80896354913783497
• MSE: 1.2 x10^7

Polynomial Fit: Using Highway-mpg as a Predictor Variable of Price.

• R-squared: 0.6741946663906514
• MSE: 2.05 x 10^7

SLR vs MLR

Usually having more variables helps the prediction, however if you do not have enough data you can run into trouble or many of the variables may just be noise and not be very useful

The MSE for MLR is Smaller than for SLR, and the $R^2$ for MLR is higher as well, MSE and $R^2$ both seem to indicate that MLR is a better fit than SLR

SLR vs Polynomial

The MSE for the Polynomial Fit is less than for the SLR and the $R^2$ is higher meaning that the Polynomial Fit is a better predictor based on 'highway-mpg' than the SLR

MLR vs Polynomial

The MSE for the MLR is smaller than for the Polynomial Fit, and the MLR also has a higher $R^2$ therefore the MLR is as better fit than the Polynomial in this case

Overall

The MLR has the lowest $R^2$ and the highest MSE, meaning that it is the best fit of the three models that have been evaluate

Model Evaluation

Model Evaluation tells us how our model works in the real wold. In-Sample evaluation does not give us an indication as to how our model performs under real lifr circumstances

Training and Test Sets

We typically split our data into a training and testing set and use to build and evaluate our model respectively

• Split into
  • 70% Training
  • 30% Testing
• Build and Train with Training Set
• Use Testing Set to evaluate model performance

sklearn gives us a function to split out data into a train and test set

from sklearn.model_selection import train_test_split

X_train, X_test, Y_train, Y_test = train_test_split(x_data,y_data,test_size=0.3,random_state=0)

Generalization Performance

Generalization Error is a measure of how well our data does at predicting previously unseen data

The error we obtain using our testing data is an approximation of this error

Cross Validation

Cross validation involves us splitting the data into folds and using a fold for testing and the remainder for training, and we make use of each combination of training, and evalution

We can use cross_val_score() to evaluate our out-of-sample evaluation

from sklearn.model_selection import cross_val_score
scores = cross_val_score(lr,X,Y,cv=n)

Where X is our predictor matrix, Y is our target, and n is our number of folds

If we want to use get the actual predicted values from our model we can do the following

from sklearn.model_selection import cross_val_predict
yhat = cross_val_predict(lr,X,Y,cv=n)

Overfitting, Underfitting, and Model Selection

The goal of model selection is to try to select the best function to fit our model, if our model is too simple we will have model that does not appropriately fit our data, whereas if we have a model that is too complex, it wil perfectly fit our testing data but will not be good at approximating new data

We need to be sure to fit the data, not the noise

It is also possible that the data we are trying to approximate cannot be fitted by polynomial at all, for example in the case of cyclic data

We can make use of the following code to look at the effect of order on our $R^2$ error for a model

rsq_test = []
order = [1,2,3,4]

for n in order:
    pr = PolynomialFeatures(degree=n)
    x_train_pr = pr.fit_transform(x_train[['feature']]
    x_test_pr = pr.fit_transform(x_test[['feature']]
    lm.fit(x_train_pr, y_train)
    rsq_test.append(lr.score(x_test_pr, y_test))

Ridge Regression

Ridge Regression prevents overfitting

Ridge regression controls the higher order parameters in our model by using a factor $\alpha$, increasing our $\alpha$ value will help us avoid overfitting until, however increasing it too far can lead to us underfitting the data

To use ridge regression we can do the following

from sklearn.linear_model import Ridge

ridge_model = Ridge(alpha=0.1)
ridge_model.fit(X,Y)
yhat=RidgeModel.predict(X)

We will usually start off with a small value of $\alpha$ such as 0.0001 and increase our value in orders of magnitude

Furthermore we can use Cross Validation to identify an optimal $\alpha$

Grid Search

Grid Search allows us to scan through multiple free parameters

Parameters such as $\alpha$ are not part of the training or fitting process. These parameters are called Hyperparameters

sklearn has a means of automatically iterating over these parameters called Grid Search

This allows us to use different hyperparameters to train our model, and select the model that provides the lowest MSE

The values of a grid search are simply a list whcih contains a dictionary for the parameters we want to modify

params = [{'param name':[value1, value2, value3, ...]}]

We can use Grid Search as follows

from sklearn.linear_model import Ridge
from sklearn.model_selection import GridSearchCV

params = [{'alpha':[0.0001, 0.01, 1, 10, 100]}, {'normalize':[True,False]}]

rr = Ridge()

grid = GridSearchCV(rr, params, cv=4)
grid.fit(X_train[['feature1','feature2',...], Y_train)

grid.best_estimator_

scores = grid.cv_results_

The scores dictionary will store the results for each hyperparameter combination given our inputs

Conclusion

For more specific examples of the different setions look at the appropriate Lab