<jupyter_start><jupyter_text>Introduction to Regression IntroductionWe will apply linear regression modules using the Python scikit-learn package:1. Understand the basics of applying regression models for prediction. 2. Evaluate the performance of regression models. 2. Apply a recipe for using scikit-learn to define, train and test machine learning models.   Overview of regression**Goal of regression**: to produce a model that represents the ‘best fit’ to some observed data.  - Typically the model is a function describing some type of curve (lines, parabolas, etc.) that is determined by a set of parameters (e.g., slope and intercept).  - “Best fit” means that there is an optimal set of parameters according to an evaluation criteria we choose.A regression models attempt to predict the value of one variable, known as the **dependent variable**, **response variable** or **label**, using the values of other variables, known as **independent variables**, **explanatory variables** or **features**. Single regression has one label used to predict one feature. Multiple regression uses two or more feature variables. In mathematical form, the **goal of regression** is: to find a function of some features $X$ which predicts the label value $y$.  This function can be written as follows:$$\hat{y} = f(X)$$The challenge in regression is to **learn** the function $f(X)$ so that the predictions of $\hat{y}$ are accurate.  We train the model to minimize the difference between our predicted $\hat{y}$ and the known label values $y$.  In fact, the entire field of **supervised learning** has this goal. Many machine learning models, including some of the latest deep learning methods, are a form of regression.  These methods often suffer from the same problems, including overfitting and mathematically unstable fitting methods.   Overview of linear regressionWe will work with linear regression models here. Linear regression are a foundational form of regression.The simplest case of linear regression is **single regression**, since there is a single feature. The function $f(X)$ is **linear in the model coefficients**. For a single vector of features $x$ the linear regression equation is written as follows:$$\hat{y} = a \cdot x + b$$The model coefficients are $a$, which we call the **slope** and $b$, which we call the **intercept**.  This is just the equation of a straight line for one variable. But, what are the best values of $a$ and $b$? In linear regression, $a$ and $b$ are chosen to minimize the squared error between the predictions and the known labels. This quantity is known as the **sum squared error** or SSE. For $n$ **training cases** the SSE is computed as follows:$$MSE = \sum_{i=1}^n \big( f(x_i) - y_i \big)^2\\= \sum_{i=1}^n \big( \hat{y}_i - y_i \big)^2\\= \sum_{i=1}^n \big( a \cdot x_i + b - y_i \big)^2$$The approach to regression that minimizes SSE is know as the **method of least squares**.  Execute a first linear regression example We will now train and evaluate a linear regression model.  In this case we will use simulated data, which means that we can compare the computed results to the known properties of the data.<jupyter_code>## Import packages
import numpy as np
import numpy.random as nr
import matplotlib.pyplot as plt
import sklearn.model_selection as ms
import sklearn.metrics as sklm
from sklearn import preprocessing
from sklearn import linear_model
import scipy.stats as ss
import seaborn as sns
import math

# make plots appear inline in the notebook
%matplotlib inline<jupyter_output><empty_output><jupyter_text>Simulate the data and plot the result. The data has the following properties:- The `x` variable is uniformly distributed between 0.0 and 10.0.- The `y` variable equals the `x` variable plus a Normally distributed random component. As a result, for the un-scalled data, the slope coefficient should be 1.0 and the intercept 0.0.<jupyter_code>nr.seed(34567)
x = np.arange(start = 0.0, stop = 10.0, step = 0.1)
y = np.add(x, nr.normal(scale = 1.0, size = x.shape[0]))

sns.regplot(x=x, y=y, fit_reg = False)
plt.xlabel('X')
plt.ylabel('Y')
plt.title('Data for regression')<jupyter_output><empty_output><jupyter_text>As expected, these data follow a straight line trend.  However, there is some dispersion of these data as a result of the addition of the Normally distributed noise.   1. Split the datasetWhen performing any type of machine learning, good data preparation is required to ensure good model performance.  The first step in preparing these data is to create **independently sampled** **training dataset** and **test data set**.  (In most cases, an independently sampled **evaluation dataset** will also be used.   In this case, no model improvement or comparison will be performed so this additional step is unnecessary.) **Information leakage** between the training and test processes: If the same data are used to train and test a machine learning model, there is a high likelihood that the model will simply be learning the training data.    In this case, the model may not **generalize** well.  A model that generalizes well produces consistent results when presented with new cases, never before encountered.  Conversely, a model with poor generalization might give unexpected results when presented with a new case. The random sub-samples of the data are created using a process called **Bernoulli sampling**.  Bernoulli sampling accepts a sample into the data subset with probability $p$.  In this case, the probability that a given case is in the training dataset is $p$.  The probability a case is in the test dataset then becomes $1-p$. The `train_test_split` function from the `sklearn.model_selection` module performed the required Bernoulli sampling.  The `train_test_split` function samples the index for the array containing the features and label values.  The code below performs this split.<jupyter_code>## Randomly sample cases to create independent training and test data
nr.seed(9988)
indx = range(len(x))
indx = ms.train_test_split(indx, test_size = 50)
x_train = np.ravel(x[indx[0]])
y_train = np.ravel(y[indx[0]])
x_test = np.ravel(x[indx[1]])
y_test = np.ravel(y[indx[1]])<jupyter_output><empty_output><jupyter_text>2. Scale numeric featuresNow that the dataset is split, the numeric feature column must be re-scaled.  The numeric range of a feature should not determine how much that feature determines the training of the machine learning model. For example, consider a data set with two features, age in years, typically measured in a few tens, and income, typically measured in tens or hundreds of thousands. There is no reason to believe that income is more important than age in some model, simply because its range of values is greater. To prevent this problem numeric features are scaled to the same range. There are many possible scaling methods. One simple method is known as **Min-Max** normalization. The data are scaled using the following formula to be in the range $\{ 0,1 \}$:$$x\_scaled_i = \frac{(x_i - Min(x))}{(Max(x) - Min(x))}$$where,$x_i $ is the ith sample value,$Min(X) $ is the minimum value of all samples,$Max(X) $ is the maximum value of all samples.In general, Min-Max normalization is a good choice for cases where the value being scaled has a complex distribution. For example, a variable with a distribution with multiple modes might be a good candidate for Min-Max normalization. Notice that the presence of a few outliers can distort the result by giving unrepresentative values of $Min(X)$ or $Max(X)$.We will use **Z-Score** normalization. Z-Score normalization transforms a variable so that it has zero mean and unit standard deviation (or variance). Z-Score normalization is performed using the following formula:$$x\_scaled_i = \frac{\big(x_i - \mu \big)}{\sigma}$$where,$\mu $ is the mean of the variable $X$,$\sigma $ is the standard deviation of the variable $X$. The code below uses the `StandardScaler` function from the `sklearn.preprocessing` package. This function computes the scaling coefficients for the training data. The resulting transformation is then applied to the training and test data using the `transform` method. The scaling transform should always be computed only on the training data, not the test or evaluation data. Generally, a numeric label does not need to be scaled. Other transformations may be required, however. Execute the code below that applies the Z-Score transformation to the training and test feature.<jupyter_code>## Scale the feature, being sure to use the scale of the training
## data not the test data. 
scaler = preprocessing.StandardScaler().fit(x_train.reshape(-1,1))
x_train = scaler.transform(x_train.reshape(-1,1))     
y_train = scaler.transform(y_train.reshape(-1,1))<jupyter_output><empty_output><jupyter_text>Train the regression modelWith the data prepared, it is time to train a regression model. This is done with the `sklearn.linear_model` package.  The steps for training most scikit-learn models are the same as used here:1. A model object is instantiated. Additional model specification can be performed at instantiation time.2. The model is fit using a numpy array of the features and the labels. In this case, there is only one feature so the `reshape` method is used to create an array of the correct dimension.Click this link to find additional [documentation on linear regression models with scikit-learn](http://scikit-learn.org/stable/modules/linear_model.htmlordinary-least-squares). Execute the code below to instantiate and fit the model.<jupyter_code>## define and fit the linear regression model
lin_mod = linear_model.LinearRegression()
lin_mod.fit(x_train.reshape(-1,1), y_train)<jupyter_output><empty_output><jupyter_text>As a first verification of this model you can print the model coefficients. These coefficients are attributes of the model object.<jupyter_code>print(lin_mod.intercept_)
print(lin_mod.coef_)<jupyter_output>[0.08516893]
[[1.03285118]]<jupyter_text>These coefficients are close to the values used in the simulation, $0.0$ and $1.0$.  Next, plot the predicted values computed from the training features.  The `predict` method is applied to the model with the training data.  A plot of the raw label values and the line of the predicted values or **scores** is then displayed.<jupyter_code>def plot_regression(x, y_score, y):
    ## Plot the result
    sns.regplot(x=x, y=y, fit_reg=False)
    plt.plot(x, y_score, c = 'red')
    plt.xlabel('X')
    plt.ylabel('Y')
    plt.title('Fit of model to test data')

y_score = lin_mod.predict(x_test.reshape(-1,1)) 

plot_regression(x_test, y_score, y_test)<jupyter_output><empty_output><jupyter_text>The red line appears to be a good fit to the data.  The errors between the scored values and the residuals appear to be minimal.  However, an objective evaluation of model performance is require.  Evaluate model performanceWith the model trained, it is time to evaluate the performance. This is done using the test dataset, so that there is no information leakage from the model training. As a first step, a set of performance metric are computed. There are many possible metrics used for the evaluation of regression models. Generally, these metrics are functions of the **residual value**, or difference between the predicted value or score and actual label value:$$r_i = f(x_i) - y_i = \hat{y}_i - y_i$$We will work with some of the more common metrics:- **Mean squared error** or MSE,  $$MSE = \frac{1}{N} \sum_{i=1}^N (f(x_i) - y_i)^2$$The mean squared error is identical to the variance of the residuals (with a slight bias). Recall that this metric is the one linear regression minimizes. Notice that mean square error is in units of the square of the label values.   - **Root mean squred error** or RMSE,  $$RMSE = \sqrt{ \frac{1}{N} \sum_{i=1}^N (f(x_i) - y_i)^2}$$The root mean squared error is identical to the standard deviation of the residuals (again, with a slight bias). Root mean square error is in the same units as the label values.  - **Mean absolute error** or MAE,$$MAE = \frac{1}{N} \sum_{i=1}^N |f(x_i) - y_i|$$  where $||$ is the absolute value operator.  The similar in interpretation to the root mean squared error. You may find this measure more intuitive since it is simply the average of the magnitude of the residuals.   - **Median absolute error**,$$Median\ Absolute\ Error = Median \big( \sum_{i=1}^N |f(x_i) - y_i| \big)$$   The median absolute error is a robust measure of the location parameter of the absolute residuals. If this measure is significantly different from the mean absolute error, it is likely that there are outliers in the residuals. - **R squared or $R^2$**, also known as the **coefficient of determination**,  $$R^2 = 1 - \frac{SS_{res}}{SS_{tot}}$$  where,   $SS_{res} = \sum_{i=1}^N r_i^2$, or the sum of the squared residuals,   $SS_{tot} = \sum_{i=1}^N y_i^2$, or the sum of the squared label values.  In other words, $R^2$ is  measure of the reduction in sum of squared values between the raw label values and the residuals. If the model has not reduced the sum of squares of the labels (a useless model!), $R^2 = 0$. On the other hand, if the model fits the data perfectly so all $r_i = 0$, then $R^2 = 1$. - **Adjusted R squared or $R^2_{adj}$** is $R^2$ adjusted for degrees of freedom in the model,$$R^2_{adj} = 1 - \frac{var(r)}{var(y)} = 1 - \frac{\frac{SS_{res}}{(n - p -1)}}{\frac{SS_{tot}}{(n-1)}}$$  where,   $var(r) = $ the variance of the residuals,   $var(y) = $ the variance of the labels,$n = $ the number of samples or cases,$p = $ number of model parameters.  The interpretation of $R^2_{adj}$ is the same as $R^2$. In many cases there will be little difference. However if the number of parameters is significant with respect to the number of cases, $R^2$ will give an overly optimistic measure of model performance. In general, the difference between $R^2_{adj}$ and $R^2$ becomes less significant as the number of cases $n$ grows. However, even for 'big data' models there can be a significant difference if there are a large number of model parameters.   ******Note:** Is it possible to get values of $R^2$ outside the range $\{ 0,1 \}$? Ordinarily no. But there are exceptions.   $R^2$ can only be greater than $1$ in degenerate cases. For example, if all label values are the same. But, in this case, you do not need a model to predict the label!What if you find your model gives an $R^2$ less than $0$? What can this possibly mean? This invariably means that there is a bug in your code and that the residuals of your model have greater dispersion than the original labels!****The code in the cell below uses functions from the `sklearn.metrics` package to compute some common metrics. There is no function for $R^2_{adj}$ in `sklearn.metrics`, but the adjustments for degrees of freedom are easily computed. Click this link to find additional [documentation on linear regression merics built into scikit-learn](http://scikit-learn.org/stable/modules/model_evaluation.htmlregression-metrics).<jupyter_code>def print_metrics(y_true, y_predicted, n_parameters):
    ## First compute R^2 and the adjusted R^2
    r2 = sklm.r2_score(y_true, y_predicted)
    r2_adj = r2 - (y_true.shape[0] - 1)/(y_true.shape[0] - n_parameters - 1) * (1 - r2)
    
    ## Print the usual metrics and the R^2 values
    print('Mean Square Error      = ' + str(sklm.mean_squared_error(y_true, y_predicted)))
    print('Root Mean Square Error = ' + str(math.sqrt(sklm.mean_squared_error(y_true, y_predicted))))
    print('Mean Absolute Error    = ' + str(sklm.mean_absolute_error(y_true, y_predicted)))
    print('Median Absolute Error  = ' + str(sklm.median_absolute_error(y_true, y_predicted)))
    print('R^2                    = ' + str(r2))
    print('Adjusted R^2           = ' + str(r2_adj))
    
print_metrics(y_test, y_score, 2)<jupyter_output>Mean Square Error      = 1.0178480188322825
Root Mean Square Error = 1.008884541873986
Mean Absolute Error    = 0.763059846639255
Median Absolute Error  = 0.6031435928079797
R^2                    = 0.8979014436983853
Adjusted R^2           = 0.7914582679796807<jupyter_text>Interpret these results:- The MSE and RMSE are as expected. The standard deviation of the simulated data is $1.0$. - The MAE and median absolute error have small values and are close together, indicating a good model fit and few significant outliers in the residuals. - The $R^2_{adj}$ and $R^2$ are both fairly close to one, indicating that the model is making useful predictions that are much better than the simple average of the label values. The residuals of a linear regression model should have an approximately Normal distribution. This condition can be easily tested using graphical methods, specifically a histogram and a Quantile-Quantile Normal plot. ******Note:** A common misconception is that the features or label of a linear regression model must have Normal distributions. This is not the case! Rather, the residuals (errors) of the model should be Normally distributed.   ****  The code below plots a kernel density plot and histogram of the residuals of the regression model.<jupyter_code>def hist_resids(y_test, y_score):
    ## first compute vector of residuals. 
    resids = np.subtract(y_test.reshape(-1,1), y_score.reshape(-1,1))
    ## now make the residual plots
    sns.displot(resids)
    plt.title('Histogram of residuals')
    plt.xlabel('Residual value')
    plt.ylabel('count')
    
hist_resids(y_test, y_score)<jupyter_output><empty_output><jupyter_text>This histogram and the kernel density plot look approximately Normal, but with some deviations. Overall, these residuals look reasonable for a real-world model.    Another useful plot is the **Quantile-Quantile Normal plot**, or Q-Q Normal plot. This plot displays quantiles of a standard Normal distribution on the horizontal axis and the quantiles of the residuals on the vertical axis. If the residuals were perfectly Normally distributed, these points would fall on a straight line. In real-world problems, you should expect the straight line relationship to be approximate.<jupyter_code>def resid_qq(y_test, y_score):
    ## first compute vector of residuals. 
    resids = np.subtract(y_test.reshape(-1,1), y_score.reshape(-1,1))
    ## now make the residual plots
    ss.probplot(resids.flatten(), plot = plt)
    plt.title('Residuals vs. predicted values')
    plt.xlabel('Quantiles of standard Normal distribution')
    plt.ylabel('Quantiles of residuals')
    
resid_qq(y_test, y_score)<jupyter_output><empty_output><jupyter_text>Notice that these points nearly fall along the straight line. This indicates that the residuals have a distribution which is approximately Normal.  We will now make one last diagnostic plot for this regression model, known as a **residual plot**.  - A plot of residuals vs. predicted values (scores) shows if there is structure in the residuals. - For an ideal regression model the variance or dispersion of the residuals should not change with the values of the predicted values. - The ideal residual plot should look like a 'fuzzy caterpillar' with no change vs. the predicted value. Any structure in this plot with change in predicted values indicates that the model fit changes with the predicted value.  - For example, if the residuals increase with predicted values the model can be said to predict only the smaller label values well.  - The opposite situation indicates that only large label values are well predicted. - Changes in the mid-range indicate that there is some nonlinear change with predicted values.  In other words, in any of these cases the model is not accurately computing the predicted values. Display and examine the residual plot for the regression model.<jupyter_code>def resid_plot(y_test, y_score):
    ## first compute vector of residuals. 
    resids = np.subtract(y_test.reshape(-1,1), y_score.reshape(-1,1))
    ## now make the residual plots
    sns.regplot(x=y_score, y=resids, fit_reg=False)
    plt.title('Residuals vs. predicted values')
    plt.xlabel('Predicted values')
    plt.ylabel('Residual')
    
resid_plot(y_test, y_score)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08165487207403375, 'fraction_numerical': 0.01286682832681744, 'mean_word_length': 4.605547850208044, 'pattern_counts': {'":': 0, '<': 43, '<?xml version=': 0, '>': 43, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '587391', 'license': 'MIT', 'n_tokens_mistral': 5686, 'n_tokens_neox': 5344, 'n_words': 2723, 'path': 'documents/intro_to_regression.ipynb', 'repo_name': 'liawbeile/liawbeile.github.io', 'size': '20207'}