<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Linear Regression
Linear regression is a method that allows us to study the statistical relationship between independent/predictor variables and dependent/response/outcome variables. It is simple to use and interpreting a linear regression model is straightforward. There are two variants
Step1:  Simple Linear Regression
Simple linear regression assumes that there is a linear relationship between two continuous variables; the predictor variable $x$ and the response variable $y$. This can be described as
Step2:  The green dots represent our artificial data i.e. the observed data. The red vertical lines indicate the errors or the residuals.  The blue line represent the regression line that best fit the observed data. This line is found by minimising the residuals. 
Multiple Linear Regression
Typically, we have more than one predictor variable (features) in our data. In such instances, we use Multiple Linear Regression
Step3:  Check data distributions.
Step4:  We can get check if we have missing values in the dataset
Step5:  No missing data. If any of our features contained missing values, we could use a heatmap to get a quick overview of where there missing data are and a rough picture on how much data is missing.
Step6:  Categorical features
Next step is to identify which features are categorical
Step7:  We see that CHAS is a Boolean feature. RAD (index of accessibility to radial highways) is also categorical since it has 9 distinct values. Let us fix them
Step8:  Since CHAS is a binary feature, we can use <a href="https
Step9:  Feature Selection
Before we can train a linear regression model, we must decide on which feature to pick as the independent variable. We must find a good predictor variable for the dependent variable. One way to find such a candidate is to create a heat map of the correlations between each of the features. Pearson's correlation can tell us as one variable increses, whether another variable increases (positive correlation), decreases (negative correlation), or stays the same (no correlation). Read more on <a href="https
Step10:  There is a strong correlation between our target feature MEDV (the house price) and RM which is the average number of rooms per house. There are also negative correlations with LSTAT and PTRATIO. This means that as one of the variables (LSTAT or PTRATIO) increases then our dependent variable (MEDV) tends to decrease and the other way around.
For now, we focus on the relationship between MEDV and RM. Let us create a linear model plot of these two attributes
Step11:  From the plot, it looks like there is a good linear fit since the error bars are not that big.
Train a Model
Before we can train a model, it is important to split the data into training set and testing set to avoid overfitting. Overfitting occurs when our model becomes too complex which may yield high accuracy on the training data but does not generalise well. 
Note that the train_test_split() method simply partitions data randomly. This means that we may get another coefficient if we run the method again. This means that the accuracy of our testing can vary depending on the random split. Therefore, a better approach would be to apply <a href="http
Step12:  It is straightforward to interpret the coefficient
Step13:  Generate Predictions
Now we use the fitted model to perform predictions on our test set.
Step14:  Let us check how far off the predicted values are from the observed values. We can visualise this via scatter plot. We can also plot a dashed line that indicates where perfect prediction values would lie.
Step15:  Create a histogram of our residuals/errors. If the residuals follow a normal distribution, then we have selected a correct model for the data set. Otherwise, we may have to reconsider the choice of linear regression model.
Step16:  Evaluate the Model
Next, we evaluate the performance of our model. The sci-learn library implements several <a href="http
<ASSISTANT_TASK:>
Python Code:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline
Explanation: Linear Regression
Linear regression is a method that allows us to study the statistical relationship between independent/predictor variables and dependent/response/outcome variables. It is simple to use and interpreting a linear regression model is straightforward. There are two variants:
Simple linear regression
One predictor variable, often denoted $x$
One response variable, often denoted $y$
Multiple linear regression
Two or more predictor variables
One dependent variable
End of explanation
from sklearn.linear_model import LinearRegression
# Create artificial data
X = np.array([   1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
Y = np.array([-1.1, 4, 1, 6, 4, 2, 8, 5, 12, 7])
# Create a model
lrm = LinearRegression()
lrm.fit(X.reshape(-1, 1), Y)
model_line = lrm.intercept_ + X * lrm.coef_[0]
fig, ax = plt.subplots(figsize=(8,4))
# Plot the data
ax.plot(X, Y, 'go') 
# Plot the regression line
ax.plot(X, model_line, 'bo-') # 
# Plot the residuals
ax.vlines(x=X, ymin=np.minimum(Y, model_line), ymax=np.maximum(Y, model_line), color='red')
Explanation: Simple Linear Regression
Simple linear regression assumes that there is a linear relationship between two continuous variables; the predictor variable $x$ and the response variable $y$. This can be described as: 
\begin{equation}
y \thickapprox \beta_0 + \beta_1 x
\end{equation}
where
$\thickapprox$ can be read as "is approximately modeled as". 
$\beta_0$ is called the intercept e.g. if $x=0$ then $y=\beta_0$
$\beta_1$ is called the slope (sometimes referred to as the coefficient). This basically describes the slope of the regression line e.g. if $x$ increases by 1 unit then $y$ increases by $\beta_1$ units.  
Together, $\beta_0$ and $\beta_1$ are called the model coefficients or parameters.
Learning the Model Coefficients
Since these coefficients are unknown to begin with, we use our data to estimate or learn them using <a href="https://www.boomer.org/c/p3/c06/c0607.html">least squares criteria</a>. Essentially, we find the <a href="https://onlinecourses.science.psu.edu/stat501/node/252">best fitting line</a> by minimising the sum of squared errors/residuals. The errors or residuals are the differences between the observed value and the estimated value. We square the difference to avoid negative values. The residual sum of squares (RSS) is defined as:
\begin{equation}
RSS = \sum_{i=1}^n{(y_i - \hat{y}_i)^2}
\end{equation}
where
$y_i$ is the observed value
$\hat{y}_i$ is the estimated/fitted value
We can illustrate this with a plot:
End of explanation
from sklearn.datasets import load_boston
boston = load_boston()
boston.keys()
print(boston['DESCR'])
df = pd.DataFrame(data=boston.data, columns=boston['feature_names'])
# Add a column for the target feature
df['MEDV'] = boston['target'] 
df.head()
df.info()
Explanation: The green dots represent our artificial data i.e. the observed data. The red vertical lines indicate the errors or the residuals.  The blue line represent the regression line that best fit the observed data. This line is found by minimising the residuals. 
Multiple Linear Regression
Typically, we have more than one predictor variable (features) in our data. In such instances, we use Multiple Linear Regression:
\begin{equation}
y \thickapprox \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \cdots + \beta_n x_n
\end{equation}
Explore the Boston House Prices Dataset
End of explanation
df.describe().transpose()
Explanation: Check data distributions.
End of explanation
df.isnull().any()
Explanation: We can get check if we have missing values in the dataset:
End of explanation
sns.heatmap(df.isnull(), cbar=False, yticklabels=False)
Explanation: No missing data. If any of our features contained missing values, we could use a heatmap to get a quick overview of where there missing data are and a rough picture on how much data is missing.
End of explanation
# Features with number of unique values
unq = { 
    column : df[column].nunique()
    for column in df.columns
}
unq
import operator
sorted(unq.items(), key=operator.itemgetter(1))
sorted(df['RAD'].unique())
Explanation: Categorical features
Next step is to identify which features are categorical:
End of explanation
df['CHAS'] = df['CHAS'].astype('category')
df['RAD'] = df['RAD'].astype('category')
df.info()
Explanation: We see that CHAS is a Boolean feature. RAD (index of accessibility to radial highways) is also categorical since it has 9 distinct values. Let us fix them:
End of explanation
from scipy import stats
stats.pointbiserialr(df['CHAS'], df['MEDV'])
Explanation: Since CHAS is a binary feature, we can use <a href="https://www.andrews.edu/~calkins/math/edrm611/edrm13.htm">Point-Biserial</a> to compute correlation between CHAS and MEDV:
End of explanation
plt.figure(figsize = (10,6))
sns.heatmap(df.corr(), annot=True, linewidths=1)
Explanation: Feature Selection
Before we can train a linear regression model, we must decide on which feature to pick as the independent variable. We must find a good predictor variable for the dependent variable. One way to find such a candidate is to create a heat map of the correlations between each of the features. Pearson's correlation can tell us as one variable increses, whether another variable increases (positive correlation), decreases (negative correlation), or stays the same (no correlation). Read more on <a href="https://stats.stackexchange.com/questions/119835/correlation-between-a-nominal-iv-and-a-continuous-dv-variable">Cross Validated</a>.
End of explanation
sns.lmplot(data=df, x='MEDV', y='RM')
Explanation: There is a strong correlation between our target feature MEDV (the house price) and RM which is the average number of rooms per house. There are also negative correlations with LSTAT and PTRATIO. This means that as one of the variables (LSTAT or PTRATIO) increases then our dependent variable (MEDV) tends to decrease and the other way around.
For now, we focus on the relationship between MEDV and RM. Let us create a linear model plot of these two attributes:
End of explanation
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
X = df['RM']
y = df['MEDV']
# Partition the data into training set and test set
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.33)
lrm = LinearRegression()
lrm.fit(X_train.values.reshape(-1, 1), y_train)
intercept = lrm.intercept_
coefficient = lrm.coef_[0]
print("Intercept:   {0}\nCoefficient: {1}".format(intercept, coefficient))
Explanation: From the plot, it looks like there is a good linear fit since the error bars are not that big.
Train a Model
Before we can train a model, it is important to split the data into training set and testing set to avoid overfitting. Overfitting occurs when our model becomes too complex which may yield high accuracy on the training data but does not generalise well. 
Note that the train_test_split() method simply partitions data randomly. This means that we may get another coefficient if we run the method again. This means that the accuracy of our testing can vary depending on the random split. Therefore, a better approach would be to apply <a href="http://scikit-learn.org/stable/modules/cross_validation.html">k-fold cross-validation</a>. However, cross-validation may not work well when the data is not uniformly distributed.
End of explanation
print('A unit increase in the RM feature is associated with a {0} units increase in the house price.'.format(coefficient))
Explanation: It is straightforward to interpret the coefficient:
End of explanation
predicted = lrm.predict(X_test.values.reshape(-1, 1))
Explanation: Generate Predictions
Now we use the fitted model to perform predictions on our test set.
End of explanation
fig, ax = plt.subplots()
ax.scatter(y_test, predicted)
ax.plot([y.min(), y.max()], [y.min(), y.max()], 'k--', linewidth=3)
ax.set_xlabel('Observed')
ax.set_ylabel('Predicted')
Explanation: Let us check how far off the predicted values are from the observed values. We can visualise this via scatter plot. We can also plot a dashed line that indicates where perfect prediction values would lie.
End of explanation
sns.distplot(y_test - predicted);
Explanation: Create a histogram of our residuals/errors. If the residuals follow a normal distribution, then we have selected a correct model for the data set. Otherwise, we may have to reconsider the choice of linear regression model.
End of explanation
from sklearn import metrics
print('Mean Absolute Error:      ', metrics.mean_absolute_error(y_test, predicted))
print('Median Absolute Error:    ', metrics.median_absolute_error(y_test, predicted))
print('Explained Variance Score: ', metrics.explained_variance_score(y_test, predicted))
print('Coefficient of Determination (R^2): ', metrics.r2_score(y_test, predicted))
print('Mean Squared Error:       ', metrics.mean_squared_error(y_test, predicted))
print('Root Mean Squared Error:', np.sqrt(metrics.mean_squared_error(y_test, predicted)))
Explanation: Evaluate the Model
Next, we evaluate the performance of our model. The sci-learn library implements several <a href="http://scikit-learn.org/stable/modules/model_evaluation.html#regression-metrics">metric functions</a>:
Mean Absolute Error is the mean of the absolute value of the errors
\begin{equation}
 \frac{1}{n} \sum_{i=1}^n \mid y_i - \hat{y}_i \mid
 \end{equation}
Median Absolute Error is interesting because it is robust to outliers. 
Explained Variance Score
\begin{equation}
 1- \frac{Var(y - \hat{y})}{Var(y)}
 \end{equation}
Coefficient of Determination ($R^2$) is often used as a tool for comparing different models. Although we want to achieve an $R^2$ value close to 1, there is no guarantee that a model with such high value will be good because $R^2$ is susceptible to overfitting.
Mean Squared Error tends to be useful with real-world data because it takes "larger" errors into account 
\begin{equation}
 \frac{1}{n} \sum_{i=1}^n (y_i - \hat{y}_i )^2
 \end{equation}
Root Mean Squared Error
End of explanation
<END_TASK>