<jupyter_start><jupyter_text>Regression Analysis- so far, we've learned the main concepts behind supervised learning and trained many models for classification tasks to predict group memberships or categorial variables- regression models are used to predict target variable on a continous scale- some important applications are:    1. predicting insurance premium    2. making weather forecast    3. predicting stock market    4. predicting housing market    5. predicting sales of a company in future month, etc.- In this notebook, we will discuss the main concepts of regression models and cover the following topics:    - Exploring and visualizing datasets    - Looking at different approaches to implement linear regression models    - Training regression models that are robust to outliers    - Evaluating regression models and diagnosing common problems    - Fitting regression models to nonlinear data     Linear regression- the goal of linear regression is to model the relationship between one or multiple features and a continuous target variable Simple linear regression- the goal of simple (**univariate**) linear regression is to model the relationship between a single feature (**explanatory variable**, *x*) and a continuous-valued **target (response variable**, *y*) - the equation of a linear model with one explanatory variable is defined as follows:    - $y = w_0 + w_1x$        - $w_0$ is the $y$ axis incercept        - $w_1$ is the weight coefficient independent variable- the goal is to learn the weights of the linear equation to describe the relationship between the independent variable and the target variable- the learned weights can be used predict responses of new samples- visusally, linear regression can be understood as finding the best-fitting straight line through the training example, as shown in the following figure![Linear Regression](./images/linear-reg-1.png)- the best-fitting line is also called the regression line- the vertical lines from the regression line to the training examples are the **offsets** or **residual** -- the errors of our prediction Multiple linear regression- generalizing the linear regression model to multiple explanatory variables    - $y = w_0 x_0 + w_1 x_1 + w_nx_x = \sum_{i=0}^{n}w_i x_i = w^T x$- visualizing 2-d, fitted hyperplane of a multiple linear regression model with two features is already difficult![Multiple linear regression](./images/multiple-linear-reg.png) - due to the challenge and limitations of visualizing multiple linear regression hyperplanes in dataset with more than 2 features, we'll focus on univariate case, using simple regression models  Exploring the Housing dataset- housing dataset contains information about houses in the suburbs of Boston collected in 1978- made freely available from UCI ML Repository or Sci-kit learn - https://github.com/scikit-learn/scikit-learn/blob/main/sklearn/datasets/data/boston_house_prices.csv- 506 examples with 14 columns- feature description:    - CRIM: Per capita crime rate by town    - ZN: Proportion of residential land zoned for lots over 25,000 sq. ft.    - INDUS: Proportion of non-retail business acres per town    - CHAS: Charles River dummy variable (= 1 if tract bounds river and 0 otherwise)    - NOX: Nitric oxide concentration (parts per 10 million)    - RM: Average number of rooms per dwelling    - AGE: Proportion of owner-occupied units built prior to 1940    - DIS: Weighted distances to five Boston employment centers    - RAD: Index of accessibility to radial highways    - TAX: Full-value property tax rate per \$10,000    - PTRATIO: Pupil-teacher ratio by town    - B: $1000(Bk – 0.63)^2$, where $Bk$ is the proportion of [people of African American descent] by town    - LSTAT: Percentage of lower status of the population    - MEDV: Median value of owner-occupied homes in $\$1000s$<jupyter_code>import pandas as pd
import numpy as np

url = 'https://raw.githubusercontent.com/scikit-learn/scikit-learn/main/sklearn/datasets/data/boston_house_prices.csv'
df = pd.read_csv(url, header=1) # column header is at row 1
df<jupyter_output><empty_output><jupyter_text>Visualize the important characteristics of a dataset- **Exploratory data analysis (EDA)** allows us to visually detect the presence of outliers, distribution of the data, and the relationships between features- let's create **scatterplot matrix** that allows us to visualize the pair-wise correlations between the different features in one place- due to space constraint, we'll use some selected columns    - feel free to explore all...<jupyter_code>import matplotlib.pyplot as plt
import seaborn as sns
cols = ['LSTAT', 'INDUS', 'NOX', 'RM', 'MEDV']
g = sns.PairGrid(df.loc[:, cols])
g.map_diag(sns.histplot)
g.map_offdiag(sns.scatterplot)
g.add_legend()<jupyter_output><empty_output><jupyter_text>- eyeball some observations:    - diagonal charts are just histogram distribution of each feature on x-axix    - there's a linear relationship between RM and MEDV, the house price    - the MEDV feature seems to be normally distributed but contains several outliers  Estimating the coefficient of a regression model- use scikit learn's LinearRegression API<jupyter_code>from sklearn.linear_model import LinearRegression
# let's just use no. of bedrooms as feature for the sake of simplicity
X = df[['RM']]
X
y = df['MEDV']
y
lr = LinearRegression()
lr.fit(X, y)
y_pred = lr.predict(X)
lr.coef_
print(f'Slope, w0: {lr.coef_[0]:.3f}')
print(f'Intercept: {lr.intercept_:.3f}')<jupyter_output>Intercept: -34.671<jupyter_text>plot the regression line- the best fitting line on the dataset<jupyter_code>sns.set_theme(color_codes=True)
fig_dims = (10, 8)
fig, ax = plt.subplots(figsize=fig_dims)
ax = sns.regplot(x="RM", y="MEDV", data=df)
ax.set(ylabel="Price in \$1000s [MEDV]", xlabel="Average number of rooms [RM]")
import locale # currency formatting
locale.setlocale(locale.LC_ALL, '') # set US locale
# let's predict the price for a 5 bedroom house
num_rooms = np.array([[5.0]])
price = lr.predict(num_rooms)
print(f'''Price of {num_rooms[0][0]:.0f} bedroom house is predicted as
          {locale.currency(price[0]*1000, grouping=True )}''')<jupyter_output>Price of 5 bedroom house is predicted as
          $10,839.92<jupyter_text>Fitting a robust regression model using RANSAC- Linear regression models can be heavily impacted by the presence of outliers- outliers can be detected and removed; requires judgement as data scientist and the domain knowledge- alternative to throwing outliers is using the RANSAC - RANdom SAmple Consensus algorithm    - fits a regression model to a subset of the data (**inliers**)- We can summarize the iterative RANSAC algorithm as follows:    1. Select a random number of examples to be inliers and fit the model.    2. Test all other data points against the fitted model and add those points that fall within a user-given tolerance to the inliers.    3. Refit the model using all inliers.    4. Estimate the error of the fitted model versus the inliers.    5. Terminate the algorithm if the performance meets a certain user-defined threshold or if a fixed number of iterations were reached; go back to step 1 otherwise.    - use RANSACRegression API of scikit-learn - https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.RANSACRegressor.html<jupyter_code>from sklearn.linear_model import RANSACRegressor
ransac = RANSACRegressor(LinearRegression(), max_trials=100,
                        min_samples=50,
                        loss='absolute_loss',
                        residual_threshold=5.0,
                        random_state=0)
ransac.fit(X, y)
# let's obtain the inliers and outliers from the fitted ransac
# plot them with the linear fit
inlier_mask = ransac.inlier_mask_
outlier_mask = np.logical_not(inlier_mask)
line_X = np.arange(3, 10, 1)
line_y_ransac = ransac.predict(line_X[:, np.newaxis])

fig_dims = (10, 8)
fig, ax = plt.subplots(figsize=fig_dims)

plt.scatter(X[inlier_mask], y[inlier_mask],
            c='steelblue', edgecolor='white',
            marker='o', label='Inliers') 
plt.scatter(X[outlier_mask], y[outlier_mask], 
            c='limegreen', edgecolor='white',
            marker='s', label='Outliers') 
plt.plot(line_X, line_y_ransac, color='black', lw=2) 
plt.xlabel('Average number of rooms [RM]') 
plt.ylabel('Price in $1000s [MEDV]') 
plt.legend(loc='upper left')
print(f'Slope: {ransac.estimator_.coef_[0]:.3f}')
print(f'Intercept: {ransac.estimator_.intercept_:.3f}')
# we get a slightly different slope and intercept compared to straight linear regression
# let's predict the price for a 5 bedroom house with ransac model
num_rooms = np.array([[5.0]])
price = ransac.predict(num_rooms)
print(f'Price of {num_rooms[0][0]:.0f} bedroom house is predicted as \
          {locale.currency(price[0]*1000, grouping=True )}')<jupyter_output>Price of 5 bedroom house is predicted as           $9,583.48<jupyter_text>Evaluating the performance of linear regression models- like supervised classifiers, regressors need to be trained on training set and evaluate on test set- the goal is to evaluate the model's performance on unseen data to estimate the generalization performance- in order to properly evaluate the model, we'll use all the variables/features in the dataset<jupyter_code>from sklearn.model_selection import train_test_split
X = df.iloc[:, :-1].values # use all the columns except for the last as explanatory variables
y = df['MEDV'].values # use the last column as dependent/response variable
# 80/20 split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=0)
slr = LinearRegression()
slr.fit(X_train, y_train)
# get the training and testing prediction
y_train_pred = slr.predict(X_train)
y_test_pred = slr.predict(X_test)<jupyter_output><empty_output><jupyter_text>Residual plots- since our model uses multiple explanatory variables, we can't visualize the linear regression line or hyperplane- residual plots lets us visualize the residual (differences or vertical distances between the actual and predicted values) versus the predicted values to diagnose our regression model- the following code plots a residual plot by simply subtracting the true target value from predicted responses<jupyter_code>fig, ax = plt.subplots(figsize=(10, 8))
plt.scatter(y_train_pred, y_train_pred - y_train, 
            c='steelblue', marker='o', 
            edgecolor='white', label='Training data')
plt.scatter(y_test_pred, y_test_pred - y_test, c='limegreen', 
            marker='s', edgecolor='white', label='Test data')
plt.xlabel('Predicted prices in $1000s')
plt.ylabel('Residuals')
plt.legend(loc='upper left')
plt.hlines(y=0, xmin=-10, xmax=50, color='black', lw=2)
plt.xlim([-10, 50])
plt.show()<jupyter_output><empty_output><jupyter_text>Mean Squared Error (MSE)- useful quantative measure of regressor models' performance- simply the averaged value of the Sum of Squared Errors (SSE)    - $MSE = \frac{1}{n}\sum_{1=1}^{n}(y^{(i)}-\hat y^{(i)})^2$<jupyter_code>from sklearn.metrics import mean_squared_error
mse_train = mean_squared_error(y_train, y_train_pred)
mse_test = mean_squared_error(y_test, y_test_pred)
print(f'MSE train {mse_train:.3f} test: {mse_test:.3f}')<jupyter_output>MSE train 19.326 test: 33.449<jupyter_text>interpreting errors- if the difference in train and test mean squared errors (MSE) is big, means the model is overfitting the training data- the interpretration of MSE depends on the dataset and feature scaling- e.g., if the house prices were represented as multiples of 1,000 (with K suffix), the same model would yield a lower MSE compared to a model that worked with unscaled features    - $(10K - 15K)^2 < (10,000 - 15,000)^2$     Coefficient of determination ($R^2$)- standarized version of the MSE for better interpretation of the MSE- $R^2$ is the fraction of response variance captured by the model          $R^2 = 1 - \frac {MSE}{Var(y)}$- higer the $R^2$, better the prediction meaning lesser the error<jupyter_code>from sklearn.metrics import r2_score
print('R^2 train: %.3f, test: %.3f' %
      (r2_score(y_train, y_train_pred), r2_score(y_test, y_test_pred)))<jupyter_output>R^2 train: 0.773, test: 0.589<jupyter_text>Running a linear regression model into a curve - polynomial regression- linearity assumption in data can be violated with polynomial regression model by adding polynomial terms:    - $y=w_0+w_1x+w_2x^2 + ... + w_nx^n$    - $n$ denotes the degree of the polynomial- although we can use ploynomial regression to model nonlinear relationshop, it is still considered a multiple linear regression model because of regression coefficients, $w$  Adding ploynomial terms- can use `PolynomialFeatures` trasformer in scikit-learn to add a quadratic term (degree = 2)- compare the linear with the polynomial fit<jupyter_code>from sklearn.preprocessing import PolynomialFeatures
# create a toy dataset
X_train = np.array([258.0, 270.0, 294.0, 320.0, 342.0,
              368.0, 396.0, 446.0, 480.0, 586.0])[:, np.newaxis]
X_train
y_train = np.array([ 236.4, 234.4, 252.8, 298.6, 314.2, 342.2, 360.8, 368.0, 391.2, 390.8])
y_train
# add a second degree polynomial term
quadratic = PolynomialFeatures(degree=2) 
X_quad = quadratic.fit_transform(X_train)
X_quad
# fit a simple linear regression model for comparison
lr = LinearRegression()
# fit linear features
lr.fit(X_train, y_train)
X_test = np.arange(250, 600, 10)[:, np.newaxis]
y_lin_test = lr.predict(X_test)
# fit quadratic features
pr = LinearRegression()
pr.fit(X_quad, y_train)
y_quad_test = pr.predict(quadratic.fit_transform(X_test))
# plot results
fig, ax = plt.subplots(figsize=(10, 5))
plt.scatter(X_train, y_train, label='Training points')
plt.plot(X_test, y_lin_test, label='Linear fit', linestyle='--')
plt.plot(X_test, y_quad_test, label='Quadratic fit')
plt.xlabel('Explanatory variable')
plt.ylabel('Predicted or known target values')
plt.legend(loc='upper left')

plt.tight_layout()
#plt.savefig('images/10_11.png', dpi=300)
plt.show()
# find the MSE and R^2
y_lin_pred = lr.predict(X_train)
y_quad_pred = pr.predict(X_quad)
print('Training MSE linear: %.3f, quadratic: %.3f' % (
        mean_squared_error(y_train, y_lin_pred),
        mean_squared_error(y_train, y_quad_pred)))
print('Training R^2 linear: %.3f, quadratic: %.3f' % (
        r2_score(y_train, y_lin_pred),
        r2_score(y_train, y_quad_pred)))<jupyter_output>Training MSE linear: 569.780, quadratic: 61.330
Training R^2 linear: 0.832, quadratic: 0.982<jupyter_text>Modeling nonlinear relationships in the Housing dataset- let's model the relationship between house prices and LSTAT (percentage of lower status of the population) using second-degree (quadratic) and third-degree (cubic) polynomials- compare quadratic and cubic polynomaials with linear fit<jupyter_code># use just one feature LSTAT as explanatory feature
X = df[['LSTAT']].values
y = df['MEDV'].values # target variable
regr = LinearRegression()
# create quadratic features
quadratic = PolynomialFeatures(degree=2)
cubic = PolynomialFeatures(degree=3)
X_quad = quadratic.fit_transform(X)
X_cubic = cubic.fit_transform(X)
# test data
X_fit = np.arange(X.min(), X.max(), 1)[:, np.newaxis]
# fit linear features
regr = regr.fit(X, y)
y_lin_fit = regr.predict(X_fit)
linear_r2 = r2_score(y, regr.predict(X))
# fit quadratic features
regr = regr.fit(X_quad, y)
y_quad_fit = regr.predict(quadratic.fit_transform(X_fit))
quadratic_r2 = r2_score(y, regr.predict(X_quad))
# fit cubic features
regr = regr.fit(X_cubic, y)
y_cubic_fit = regr.predict(cubic.fit_transform(X_fit))
cubic_r2 = r2_score(y, regr.predict(X_cubic))
# plot results
fig, ax = plt.subplots(figsize=(10, 5))
plt.scatter(X, y, label='Training points', color='lightgray')

plt.plot(X_fit, y_lin_fit, 
         label='Linear (d=1), $R^2=%.2f$' % linear_r2, 
         color='blue', 
         lw=2, 
         linestyle=':')

plt.plot(X_fit, y_quad_fit, 
         label='Quadratic (d=2), $R^2=%.2f$' % quadratic_r2,
         color='red', 
         lw=2,
         linestyle='-')

plt.plot(X_fit, y_cubic_fit, 
         label='Cubic (d=3), $R^2=%.2f$' % cubic_r2,
         color='green', 
         lw=2, 
         linestyle='--')

plt.xlabel('% lower status of the population [LSTAT]')
plt.ylabel('Price in $1000s [MEDV]')
plt.legend(loc='upper right')

plt.show()<jupyter_output><empty_output><jupyter_text>other transformations ?- polynomial is not always the best choice for modeling non-linear relationships- e.g., MEDV-LSTAT scatterplot may lead to a hypothesis that a log-transformation of the LSTAT feature variable and the square root of MEDV may project the data onto linear feature space suitable for a linear regression fit    - $f(x) = e^{-x}$    - $log(f(x)) = -x$- natural log of an exponential function is a straight line<jupyter_code>X = df[['LSTAT']].values
y = df['MEDV'].values

# transform features
X_log = np.log(X)
y_sqrt = np.sqrt(y)

# fit features
X_fit = np.arange(X_log.min()-1, X_log.max()+1, 1)[:, np.newaxis]

regr = regr.fit(X_log, y_sqrt)
y_lin_fit = regr.predict(X_fit)
linear_r2 = r2_score(y_sqrt, regr.predict(X_log))

# plot results
plt.scatter(X_log, y_sqrt, label='Training points', color='lightgray')

plt.plot(X_fit, y_lin_fit, 
         label='Linear (d=1), $R^2=%.2f$' % linear_r2, 
         color='blue', 
         lw=2)

plt.xlabel('log(% lower status of the population [LSTAT])')
plt.ylabel('$\sqrt{Price \; in \; \$1000s \; [MEDV]}$')
plt.legend(loc='lower left')

plt.tight_layout()
#plt.savefig('images/10_13.png', dpi=300)
plt.show()<jupyter_output><empty_output><jupyter_text>Decision Tree Regressor- Decision Tree can be used as a regressor model![DT Regressor](./images/DT-Regressor.png)- points to note:    1. no data transformation required; feature is analyzed one at a time    2. DT regressor captures the general trend in the data    3. doesn't capture the contuinuity and differentiability of the desired predition    4. need to be careful about choosing an appropriate value for the depth of the tree so as to not overfit or underfit the data    - let's visualize fitting regression curve with decision tree<jupyter_code>from sklearn.tree import DecisionTreeRegressor
def lin_regplot(X, y, model):
    plt.scatter(X, y, c='steelblue', edgecolor='white', s=70)
    plt.plot(X, model.predict(X), color='black', lw=2)    
    return 
X = df[['LSTAT']].values
y = df.MEDV
tree = DecisionTreeRegressor(max_depth=3)
tree.fit(X, y)
sort_idx = X.flatten().argsort()
fig, ax = plt.subplots(figsize=(10, 5))
lin_regplot(X[sort_idx], y[sort_idx], tree)
plt.xlabel('% lower status of the population [LSTAT]')
plt.ylabel('Price in $1000s [MEDV]')
plt.show()
# making predictions
# generate 10 random LSTAT samples using np.random
X_test = np.random.choice(X.max(axis=1), 10) # random values between 0 and max LSTAT
X_test = X_test[:, np.newaxis]
X_test
y_pred = tree.predict(X_test)
y_pred<jupyter_output><empty_output><jupyter_text>MSE/$R^2$? - since we don't know the actual y for the random test data, we can't calculate the MSE/$R^2$- let's split the dataset into train/test and evaluate the DT Regressor model<jupyter_code># 80/20 split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=0)
tree = DecisionTreeRegressor(max_depth=3)
tree.fit(X_train, y_train)
y_pred = tree.predict(X_test)
y_train_pred = tree.predict(X_train)
print('Training MSE linear: %.3f'%
        mean_squared_error(y_train, y_train_pred))
print('Training R^2 linear: %.3f'%
        r2_score(y_train, y_train_pred))
print('Testing MSE linear: %.3f'%
        mean_squared_error(y_test, y_pred))
print('Testing R^2 linear: %.3f'%
        r2_score(y_test, y_pred))<jupyter_output>Testing MSE linear: 39.417
Testing R^2 linear: 0.516<jupyter_text>Underfitting and Overfitting Decision Tree- DT Regressor models can suffer from underfit and overfit- one of the reasons could be "Tree Depth"- higher the the tree depth, more leaf nodes    - tree with 10 levels, will have $2^{10} = 1024$ leaves- leaves with very fewer houses will make predictions that are quite close to those homes' actual values    - however, they may make unreliable predictions for new data    - because, each prediction is based on only a few houses    ![DT under and overfitting](./images/dt-fitting-errors.png)  Finding the optimal leaf nodes- Scit-kit learn DT Regressor provides `max_leaf_nodes` argument to control overfitting vs underfitting- we can use a utility function to help compare MSE scores from models with different `max_leaf_nodes`<jupyter_code>def get_mse(max_leaf_nodes, X_train, y_train, X_test, y_test):
    model = DecisionTreeRegressor(max_leaf_nodes=max_leaf_nodes, random_state=0)
    model.fit(X_train, y_train)
    y_preds = model.predict(X_test)
    mse = mean_squared_error(y_test, y_preds)
    return mse
# compare MAE with differing values of max_leaf_nodes
errors = []
for i in range(1, 11, 1):
    max_leaf_nodes = 2**i
    mse_err = get_mse(max_leaf_nodes, X_train, y_train, X_test, y_test)
    print("Max leaf nodes: %d  \t\t Mean Squared Error:  %d" %(max_leaf_nodes, mse_err))
    errors.append((mse_err, max_leaf_nodes))
mse, optimal_leaves = min(errors)
print(f'Optimal leaves: {optimal_leaves} MSE: {mse:.3f}')
# let's find the optimal depth of the tree
import math
print(f'Optimal depth = {math.floor(math.log2(optimal_leaves))}')<jupyter_output>Optimal depth = 3<jupyter_text>Random forest regression- a random forest is an ensemble technique that combines multiple decision trees- a random forest usually has a better generalization performance than an individual decision tree due to randomness (which helps to decrease the model's variance)- RF are also less sensitive to outliers and don't require much parameter tuning- only parameter that could be tuned is number of trees in the forest- https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.RandomForestRegressor.html- let's apply `RandomForestRegression` to Housing dataset<jupyter_code># use all the features
X = df.iloc[:, :-1].values
y = df.MEDV
X_train, X_test, y_train, y_test = train_test_split(X, y, 
                                    test_size=.2, random_state=1)
X_train, X_test, y_train, y_test = train_test_split(X, y, 
                                    test_size=.2, random_state=1)
from sklearn.ensemble import RandomForestRegressor
forest = RandomForestRegressor(random_state=1, n_jobs=-1)
# train and test the model and find MSE and R^2
forest.fit(X_train, y_train)
y_train_pred = forest.predict(X_train)
y_test_pred = forest.predict(X_test)
print('MSE train: %.3f, test: %.3f' % (
        mean_squared_error(y_train, y_train_pred),
        mean_squared_error(y_test, y_test_pred)))
print('R^2 train: %.3f, test: %.3f' % (
        r2_score(y_train, y_train_pred),
        r2_score(y_test, y_test_pred)))
# Let's plot the residuals of the predictions
fig, ax = plt.subplots(figsize=(10, 6))
plt.scatter(y_train_pred,  
            y_train_pred - y_train, 
            c='steelblue',
            edgecolor='white',
            marker='o', 
            s=35,
            alpha=0.9,
            label='Training data')
plt.scatter(y_test_pred,  
            y_test_pred - y_test, 
            c='limegreen',
            edgecolor='white',
            marker='s', 
            s=35,
            alpha=0.9,
            label='Test data')

plt.xlabel('Predicted prices in $1000s')
plt.ylabel('Residual values in $1000s')
plt.legend(loc='upper left')
plt.hlines(y=0, xmin=-10, xmax=50, lw=2, color='black')
plt.xlim([-10, 50])
plt.tight_layout()

plt.show()<jupyter_output><empty_output><jupyter_text>Feature ranking- both `RandomForestClassifier` and `RandomForestRegressor` provide `feature_importances_` attribute- the following code displays the ranking and bar charts of each feature based on its importance value- model trained above is used to demostrate it<jupyter_code>feat_labels = df.columns[:-1]
importances = forest.feature_importances_

# https://numpy.org/doc/stable/reference/generated/numpy.argsort.html
# return the indices that would sort the importances array and reverse it
indices = np.argsort(importances)[::-1]

# print all the features and their importances in highest to lowest importance
for f in range(X_train.shape[1]):
    print("%2d) %-*s %f" % (f + 1, 30, 
                            feat_labels[indices[f]], 
                            importances[indices[f]]))

# plot the histogram bar chart
plt.title('Housing Dataset - Feature Importance')
plt.bar(range(X_train.shape[1]), 
        importances[indices],
        align='center')

plt.xticks(range(X_train.shape[1]), 
           feat_labels[indices], rotation=90)
plt.xlim([-1, X_train.shape[1]])
plt.tight_layout()
#plt.savefig('images/04_09.png', dpi=300)
plt.show()
indices<jupyter_output><empty_output><jupyter_text>Use top features to train RandomForestRegressor- use forward selection technique from highest to lowest important feature- can also use to train and test other regressor models<jupyter_code>df
MSEs = [] # collect all the MSE
R2s = [] # collect all the R^2
feature_ids = []
y = df.MEDV
for feature_id in indices:
    feature_ids.append(feature_id)
    X = df.iloc[:, feature_ids].values
    X_train, X_test, y_train, y_test = train_test_split(X, y, 
                                    test_size=.2, random_state=1)
    forest = RandomForestRegressor(random_state=1, n_jobs=-1)
    forest.fit(X_train, y_train)
    y_train_pred = forest.predict(X_train)
    y_test_pred = forest.predict(X_test)
    mse = mean_squared_error(y_test, y_test_pred)
    print('Top features: %d MSE train: %.3f, test: %.3f' % (
        len(feature_ids),
        mean_squared_error(y_train, y_train_pred),
        mse))
    r2 = r2_score(y_test, y_test_pred)
    print('Top features: %d R^2 train: %.3f, test: %.3f' % (
        len(feature_ids),
        r2_score(y_train, y_train_pred),
        r2))
    MSEs.append((mse, len(feature_ids)))
    R2s.append((r2, len(feature_ids)))
# let's see the MSEs for test datasets
MSEs
# can plot 
import matplotlib.pyplot as plt
data = pd.DataFrame(MSEs)
plt.plot(data[1], data[0])
plt.xlabel("Top Feature Count")
plt.ylabel("Mean Squared Error (MSE)")
plt.title("Feature ranking on Housing Dataset")
plt.show()
# let's see the R2s for test datasets
R2s
# can do line plot as well<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup