# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# + [markdown] id="omJl3NRs12zw" colab_type="text"
# _Lambda School Data Science_
#
# # Ordinary Least Squares Regression - notes
#
# ## What is Linear Regression?
#
# Linear Regression is a statistical model that seeks to describe the relationship between some y variable and one or more x variables. 
#
# ![Linear Regression](https://upload.wikimedia.org/wikipedia/commons/thumb/3/3a/Linear_regression.svg/1200px-Linear_regression.svg.png)
#
# In the simplest case, linear regression seeks to fit a straight line through a cloud of points. This line is referred to as the "regression line" or "line of best fit." This line tries to summarize the relationship between our X and Y in a way that enables us to use the equation for that line to make predictions.
#
# ### Synonyms for "y variable" 
# - Dependent Variable
# - Response Variable
# - Outcome Variable 
# - Predicted Variable
# - Measured Variable
# - Explained Variable
# - Label
# - Target
#
# ### Synonyms for "x variable"
# - Independent Variable
# - Explanatory Variable
# - Regressor
# - Covariate
# - Feature

# + [markdown] id="GLfiyLr-12z0" colab_type="text"
# # Simple Linear Regresion (bivariate)
#
# ## Making Predictions
#
# Say that we were trying to create a model that captured the relationship between temperature outside and ice cream sales. In Machine Learning our goal is often different that of other flavors of Linear Regression Analysis, because we're trying to fit a model to this data with the intention of making **predictions** on new data (in the future) that we don't have yet.
#
# ## What are we trying to predict?
#
# So if we had measured ice cream sales and the temprature outside on 11 different days, at the end of our modeling **what would be the thing that we would want to predict? - Ice Cream Sales or Temperature?**
#
# We would probably want to be measuring temperature with the intention of using that to **forecast** ice cream sales. If we were able to successfully forecast ice cream sales from temperature, this might help us know beforehand how much ice cream to make or how many cones to buy or on which days to open our store, etc. Being able to make predictions accurately has a lot of business implications. This is why making accurate predictions is so valuable (And in large part is why data scientists are paid so well).
#
# ### Y Variable Intuition
#
# We want the thing that we're trying to predict to serve as our **y** variable. This is why it's sometimes called the "predicted variable." We call it the "dependent" variable because our prediction for how much ice cream we're going to sell "depends" on the temperature outside. 
#
# ### X Variable Intuition
#
# All other variables that we use to predict our y variable (we're going to start off just using one) we call our **x** variables. These are called our "independent" variables because they don't *depend* on y, they "explain" y. Hence they are also referred to as our "explanatory" variables.

# + id="p0X_but-12z5" colab_type="code" colab={}
# %matplotlib inline
from ipywidgets import interact
from matplotlib.patches import Rectangle
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score
import statsmodels.api as sm

# + id="i17Ht2zf120L" colab_type="code" colab={}
columns = ['Year','Incumbent Party Candidate','Other Candidate','Incumbent Party Vote Share']

data = [[1952,"Stevenson","Eisenhower",44.6],
        [1956,"Eisenhower","Stevenson",57.76],
        [1960,"Nixon","Kennedy",49.91],
        [1964,"Johnson","Goldwater",61.34],
        [1968,"Humphrey","Nixon",49.60],
        [1972,"Nixon","McGovern",61.79],
        [1976,"Ford","Carter",48.95],
        [1980,"Carter","Reagan",44.70],
        [1984,"Reagan","Mondale",59.17],
        [1988,"<NAME>r.","Dukakis",53.94],
        [1992,"<NAME>r.","Clinton",46.55],
        [1996,"Clinton","Dole",54.74],
        [2000,"Gore","Bush, Jr.",50.27],
        [2004,"Bush, Jr.","Kerry",51.24],
        [2008,"McCain","Obama",46.32],
        [2012,"Obama","Romney",52.00]]
        
df = pd.DataFrame(data=data, columns=columns)

# + id="TmAj9NrL120X" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 545} outputId="6851ba3d-9ecc-46f2-daf6-8884ade0a64a"
df

# + id="pT0AAnic120j" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 283} outputId="330a08e4-145a-4cf9-b245-3ee08d8f52a1"
df.plot(x='Year', y='Incumbent Party Vote Share');

# + id="Y8MKji6k120z" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 170} outputId="3cae751c-3d9d-4e0e-efec-a1c799ef4745"
df['Incumbent Party Vote Share'].describe()

# + id="IuMStnkH1203" colab_type="code" colab={}
target = 'Incumbent Party Vote Share'
df['Prediction'] = df[target].mean()
df['Error'] = df['Prediction'] - df[target]

# + id="dBJlc_v3120_" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 545} outputId="a7e2c6c8-e451-4ca2-bdda-4c5330ffce27"
df

# + id="AhqQBZsV121F" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 34} outputId="85704f76-bb97-43b4-eeb7-e69115b9680c"
df['Error'].sum()

# + id="GgbrVyC2121N" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 34} outputId="0387cacd-8081-4ddb-b8da-f6a5215af385"
df['Absolute Error'] = df['Error'].abs()
df['Absolute Error'].sum()

# + id="wPmTNGtk121S" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 34} outputId="a7516d03-61d2-4454-cd9b-b07e2765f3e8"
df['Absolute Error'].mean()

# + id="00PO9ett121W" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 34} outputId="4ccc2774-1327-4720-9eb0-898eb9e8b42a"
mean_absolute_error(y_true=df[target], y_pred=df['Prediction'])

# + id="YaSlOKCZ121b" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 34} outputId="dde2d83f-9cfb-4595-943e-a9dd15a66b91"
r2_score(y_true=df[target], y_pred=df['Prediction'])

# + id="hcQLXMEF121h" colab_type="code" colab={}
columns = ['Year','Average Recent Growth in Personal Incomes']

data = [[1952,2.40],
        [1956,2.89],
        [1960, .85],
        [1964,4.21],
        [1968,3.02],
        [1972,3.62],
        [1976,1.08],
        [1980,-.39],
        [1984,3.86],
        [1988,2.27],
        [1992, .38],
        [1996,1.04],
        [2000,2.36],
        [2004,1.72],
        [2008, .10],
        [2012, .95]]
        
growth = pd.DataFrame(data=data, columns=columns)

# + id="MWe0Nj9A121l" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 562} outputId="8bb4f5d9-ca2b-43b0-dbc9-e0a044bbda17"
df = df.merge(growth)
df

# + id="h-Cyego5121p" colab_type="code" colab={}
feature = 'Average Recent Growth in Personal Incomes'

# + id="L4jA1Grm121u" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 283} outputId="560a6e82-fe0b-4957-ba11-e1191981482d"
df.plot(x=feature, y=target, kind='scatter');

# + [markdown] id="u7OoSyAI121y" colab_type="text"
# We can see from the scatterplot that these data points seem to follow a somewhat linear relationship. This means that we could probably summarize their relationship well by fitting a line of best fit to these points. Lets do it.
#
#
# ## The Equation for a Line
#
# As we know a common equation for a line is:
#
# \begin{align}
# y = mx + b
# \end{align}
#
# Where $m$ is the slope of our line and $b$ is the y-intercept. 
#
# If we want to plot a line through our cloud of points we figure out what these two values should be. Linear Regression seeks to **estimate** the slope and intercept values that describe a line that best fits the data points.

# + id="OuoILQfr1210" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 34} outputId="e17acaa4-fff5-45a8-f52d-b78cc402693d"
m = 4
b = 44
df['Prediction'] = m * df[feature] + b
df['Error'] = df['Prediction'] - df[target]
df['Absolute Error'] = df['Error'].abs()
df['Absolute Error'].sum()

# + id="0_k2PYLD1215" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 34} outputId="924b038f-914a-4691-fdcb-e0a1e741856c"
df['Absolute Error'].mean()

# + id="J73RTrlJ1218" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 34} outputId="0345544b-2090-4041-8f9c-71416f9d8a8b"
r2_score(y_true=df[target], y_pred=df['Prediction'])

# + id="-LVNdAdL122A" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 283} outputId="3c5bf43f-5f8d-47a5-cc29-28dc8ce951b1"
ax = df.plot(x=feature, y=target, kind='scatter')
df.plot(x=feature, y='Prediction', kind='line', ax=ax);


# + id="RkyVu4qE122H" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 458} outputId="1a5f0003-5571-4a8f-8287-3b12d27efe1c"
def regression(m, b):
    df['Prediction'] = m * df[feature] + b
    df['Error'] = df['Prediction'] - df[target]
    df['Absolute Error'] = df['Error'].abs()
    sum_absolute_error = df['Absolute Error'].sum()
    
    title = f'Sum of absolute errors: {sum_absolute_error}'
    ax = df.plot(x=feature, y=target, kind='scatter', title=title, figsize=(7, 7))
    df.plot(x=feature, y='Prediction', kind='line', ax=ax)
    
regression(m=3, b=46)


# + [markdown] id="fl7lVFED122L" colab_type="text"
# ## Residual Error 
#
# The residual error is the distance between points in our dataset and our regression line.

# + id="x3F3aMHc122M" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 458} outputId="dc601b31-98c6-481e-9c3f-4cba4319415d"
def regression(m, b):
    df['Prediction'] = m * df[feature] + b
    df['Error'] = df['Prediction'] - df[target]
    df['Absolute Error'] = df['Error'].abs()
    sum_absolute_error = df['Absolute Error'].sum()
    
    title = f'Sum of absolute errors: {sum_absolute_error}'
    ax = df.plot(x=feature, y=target, kind='scatter', title=title, figsize=(7, 7))
    df.plot(x=feature, y='Prediction', kind='line', ax=ax)
    
    for x, y1, y2 in zip(df[feature], df[target], df['Prediction']):
        ax.plot((x, x), (y1, y2), color='grey')

regression(m=3.5, b=46)

# + id="XoYM68Zu122P" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 17} outputId="5a421d96-2c14-4120-c6df-e3fbfa3b0515"
interact(regression, m=(-10,10,0.5), b=(40,60,0.5));

# + id="4QPg_F2U122R" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 562} outputId="2a6c51e0-73ac-4417-908d-acc99fc4101e"
df

# + id="yLcIg0Df122T" colab_type="code" colab={}
df['Square Error'] = df['Error'] **2

# + id="vyc477ez122X" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 562} outputId="5649e29e-e014-46b0-d1f1-a27949b6fa65"
df


# + id="h5sDSbXI122a" colab_type="code" colab={}
def regression(m, b):
    df['Prediction'] = m * df[feature] + b
    df['Error'] = df['Prediction'] - df[target]
    df['Absolute Error'] = df['Error'].abs()
    df['Square Error'] = df['Error'] **2
    sum_square_error = df['Square Error'].sum()
    
    title = f'Sum of square errors: {sum_square_error}'
    ax = df.plot(x=feature, y=target, kind='scatter', title=title, figsize=(7, 7))
    df.plot(x=feature, y='Prediction', kind='line', ax=ax)
        
    xmin, xmax = ax.get_xlim()
    ymin, ymax = ax.get_ylim()
    scale = (xmax-xmin)/(ymax-ymin)

    for x, y1, y2 in zip(df[feature], df[target], df['Prediction']):
        bottom_left = (x, min(y1, y2))
        height = abs(y1 - y2)
        width = height * scale
        ax.add_patch(Rectangle(xy=bottom_left, width=width, height=height, alpha=0.1))


# + id="KiCmd3hP122d" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 17} outputId="9322cc3b-d3a8-403f-edd4-72429e0d7c1a"
interact(regression, m=(-10,10,0.5), b=(40,60,0.5));

# + id="-WZCy-mkAV7V" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 458} outputId="e608b88b-f721-4274-c357-4137d4bf491b"
regression(3, 46.5)

# + id="kyR-P28w122g" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 1289} outputId="62f84c38-ebbc-4e8c-c9c2-148c9550df4c"
b = 46
ms = np.arange(-10,10,0.5)
sses = []

for m in ms:
    predictions = m * df[feature] + b
    errors = predictions - df[target]
    square_errors = errors ** 2
    sse = square_errors.sum()
    sses.append(sse)
    
hypotheses = pd.DataFrame({'Slope': ms})
hypotheses['Intercept'] = b
hypotheses['Sum of Square Errors'] = sses

hypotheses

# + id="2AnrMiRP122j" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 295} outputId="392e6926-d65a-439e-873b-01c8fc06312d"
hypotheses.plot(x='Slope', y='Sum of Square Errors', 
                title=f'Intercept={b}');

# + id="0MXhGpft122o" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 34} outputId="168fc54f-a239-47cf-c5a4-b5258ca8bde4"
X = df[[feature]]
y = df[target]

X.shape, y.shape

# + id="nVPrfuRt122r" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 51} outputId="b577bf80-94f1-4f5d-bca2-4be41a0e2970"
model = LinearRegression()
model.fit(X, y)

# + id="8oo33tpb122v" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 34} outputId="79ba6bd3-f09b-42f1-eb4d-baf43c348053"
model.coef_, model.intercept_

# + id="k-tOWY_L122y" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 34} outputId="823cdbff-cc93-4ec2-e7bc-5c7db9b69222"
model.predict([[1]])

# + id="xTmPs7BO1222" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 34} outputId="6a0cab91-1153-4fcf-e25c-8d35fa0cd353"
model.predict([[2]])

# + id="tsAlmUAK1224" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 34} outputId="6841f462-ef3f-4d8f-c453-edca1cd5f68d"
model.predict([[3]])

# + id="GYGftGc81226" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 85} outputId="4d4222bc-d70a-4b8c-9a24-3fc1add57b9a"
model.predict(X)

# + id="nlGPNhKU1227" colab_type="code" colab={}
df['Prediction'] = model.predict(X)

# + id="LkMzv9SN122-" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 295} outputId="19aaabc2-22c8-4755-c3e5-187c654620ec"
ax = df.plot(x=feature, y=target, kind='scatter', title='sklearn LinearRegression')
df.plot(x=feature, y='Prediction', kind='line', ax=ax);

# + id="ApuX9mlU123C" colab_type="code" colab={}
df['Error'] = df['Prediction'] - y

# + id="Xa-zVZl_123E" colab_type="code" colab={}
df['Absolute Error'] = df['Error'].abs()
df['Square Error'] = df['Error'] ** 2

# + id="9X1kA7Ri123G" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 34} outputId="4b37b51b-92dd-45c5-e222-c5c0cad961c4"
df['Square Error'].mean()

# + id="QGIFHfZj123I" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 34} outputId="1fb992bf-d6b6-4835-a655-edea5eae2756"
mean_squared_error(y_true=y, y_pred=model.predict(X))

# + id="5nzJ8wf-123L" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 34} outputId="a87c6043-43e2-4d08-97c4-5294136206ed"
np.sqrt(mean_squared_error(y, model.predict(X)))

# + [markdown] id="W9MRIt4c123O" colab_type="text"
# # R Squared:  $R^2$
#
# One final attribute of linear regressions that we're going to talk about today is a measure of goodness of fit known as $R^2$ or R-squared. $R^2$ is a statistical measure of how close the data are fitted to our regression line. A helpful interpretation for the $R^2$ is the percentage of the dependent variable that is explained by the model.
#
# In other words, the $R^2$ is the percentage of y that is explained by the x variables included in the model. For this reason the $R^2$ is also known as the "coefficient of determination," because it explains how much of y is explained (or determined) by our x varaibles. We won't go into the calculation of $R^2$ today, just know that a higher $R^2$ percentage is nearly always better and indicates a model that fits the data more closely. 

# + id="dSDVxwF3123P" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 34} outputId="f17fd758-bc00-4f2f-a208-4b22b6b1b785"
model.score(X, y)

# + id="osXOE4Ew123R" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 34} outputId="4b145b79-485e-4296-e359-5ae87244db73"
r2_score(y, model.predict(X))

# + [markdown] id="Rsrdnk1d123V" colab_type="text"
# ### Statsmodels

# + [markdown] id="0zt4dl8O123V" colab_type="text"
# https://www.statsmodels.org/dev/examples/notebooks/generated/ols.html

# + id="_Jov649z123W" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 550} outputId="1fb1986b-778b-41cf-a571-2ea14b2d7e80"
model = sm.OLS(y, sm.add_constant(X))
results = model.fit()
print(results.summary())

# + [markdown] id="8_vwwnCg123X" colab_type="text"
# # The Anatomy of Linear Regression
#
# - Intercept: The $b$ value in our line equation $y=mx+b$
# - Slope: The $m$ value in our line equation $y=mx+b$. These two values together define our regression line.
#
# ![Slope and Intercept](http://www.ryanleeallred.com/wp-content/uploads/2018/08/linear-regression-diagram.png)
#
# - $\hat{y}$ : A prediction
# - Line of Best Fit (Regression Line)
# - Predicted (fitted) Values: Points on our regression line
# - Observed Values: Points from our dataset
# - Error: The distance between predicted and observed values.
#
# ![Residual Error](http://www.ryanleeallred.com/wp-content/uploads/2018/08/residual-or-error.gif)
#

# + [markdown] id="5aYmOfe_123Z" colab_type="text"
# # More Formal Notation
#
# ![Simple Linear Regression](http://www.ryanleeallred.com/wp-content/uploads/2018/08/simple-regression-formula.png)
#
# We have talked about a line of regression being represented like a regular line $y=mx+b$ but as we get to more complicated versions we're going to need to extend this equation. So lets establish the proper terminology.
#
# **X** - Independent Variable, predictor variable, explanatory variable, regressor, covariate
#
# **Y** - Response variable, predicted variable, measured vairable, explained variable, outcome variable
#
# $\beta_0$ - "Beta Naught" or "Beta Zero", the intercept value. This is how much of y would exist if X were zero. This is sometimes represented by the letter "a" but I hate that. So it's "Beta 0" during my lecture.
#
# $\beta_1$ - "Beta One" The primary coefficient of interest. This values is the slope of the line that is estimated by "minimizing the sum of the squared errors/residuals" - We'll get to that. 
#
# $\epsilon$ - "Epsilon" The "error term", random noise, things outside of our model that affect y.

# + [markdown] id="tNz-iXN0123a" colab_type="text"
# # How Does it do it?
#
# ## Minimizing the Sum of the Squared Error
#
# The most common method of estimating our $\beta$ parameters  is what's known as "Ordinary Least Squares" (OLS). (There are different methods of arriving at a line of best fit). OLS estimates the parameters that minimize the squared distance between each point in our dataset and our line of best fit. 
#
# \begin{align}
# SSE = \sum(y_i - \hat{y})^2
# \end{align}
#
#

# + [markdown] id="lGzEoymn123a" colab_type="text"
# ## Linear Algebra!
#
# The same result that is found by minimizing the sum of the squared errors can be also found through a linear algebra process known as the "Least Squares Solution:"
#
# ![OLS Regression](http://www.ryanleeallred.com/wp-content/uploads/2018/08/OLS-linear-algebra.png)
#
# Before we can work with this equation in its linear algebra form we have to understand how to set up the matrices that are involved in this equation. 

# + [markdown] id="ayDaqc0D123b" colab_type="text"
# ### The $\beta$ vector
#
# The $\beta$ vector represents all the parameters that we are trying to estimate, our $y$ vector and $X$ matrix values are full of data from our dataset. The $\beta$ vector holds the variables that we are solving for: $\beta_0$ and $\beta_1$
#
# Now that we have all of the necessary parts we can set them up in the following equation:
#
# \begin{align}
# y = X \beta + \epsilon
# \end{align}
#
# Since our $\epsilon$ value represents **random** error we can assume that it will equal zero on average.
#
# \begin{align}
# y = X \beta
# \end{align}
#
# The objective now is to isolate the $\beta$ matrix. We can do this by pre-multiplying both sides by "X transpose" $X^{T}$.
#
# \begin{align}
# X^{T}y =  X^{T}X \beta
# \end{align}
#
# Since anything times its transpose will result in a square matrix, if that matrix is then an invertible matrix, then we should be able to multiply both sides by its inverse to remove it from the right hand side. (We'll talk tomorrow about situations that could lead to $X^{T}X$ not being invertible.)
#
# \begin{align}
# (X^{T}X)^{-1}X^{T}y =  (X^{T}X)^{-1}X^{T}X \beta
# \end{align}
#
# Since any matrix multiplied by its inverse results in the identity matrix, and anything multiplied by the identity matrix is itself, we are left with only $\beta$ on the right hand side:
#
# \begin{align}
# (X^{T}X)^{-1}X^{T}y = \hat{\beta}
# \end{align}
#
# We will now call it "beta hat" $\hat{\beta}$ because it now represents our estimated values for $\beta_0$ and $\beta_1$
#
# ### Lets calculate our $\beta$ coefficients with numpy!

# + id="6lS8o1Tx123b" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 958} outputId="661db205-ded2-4287-cf78-654a5d3fc58b"
X = sm.add_constant(df[feature]).values
print('X')
print(X)

y = df[target].values[:, np.newaxis]
print('y')
print(y)

X_transpose = X.T
print('X Transpose')
print(X_transpose)

X_transpose_X = X_transpose @ X
print('X Transpose X')
print(X_transpose_X)

X_transpose_X_inverse = np.linalg.inv(X_transpose_X)
print('X Transpose X Inverse')
print(X_transpose_X_inverse)

X_transpose_y = X_transpose @ y
print('X Transpose y')
print(X_transpose_y)

beta_hat = X_transpose_X_inverse @ X_transpose_y
print('Beta Hat')
print(beta_hat)

# + [markdown] id="eXlj3Juu123d" colab_type="text"
# # Multiple Regression
#
# Simple or bivariate linear regression involves a single $x$ variable and a single $y$ variable. However, we can have many $x$ variables. A linear regression model that involves multiple x variables is known as **Multiple** Regression - NOT MULTIVARIATE!
#
# ![Multiple Regression](http://www.ryanleeallred.com/wp-content/uploads/2018/08/multiple-regression-model.png)

# + id="yLBSiVgM123e" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 579} outputId="7e8d4079-296a-4392-b4ba-7f391e2a0790"
df.sort_values(by='Square Error', ascending=False)

# + id="YG6lPV8b123f" colab_type="code" colab={}
"""
Fatalities denotes the cumulative number of American military
fatalities per millions of US population the in Korea, Vietnam,
Iraq and Afghanistan wars during the presidential terms
preceding the 1952, 1964, 1968, 1976 and 2004, 2008 and
2012 elections.

http://www.douglas-hibbs.com/HibbsArticles/HIBBS-PRESVOTE-SLIDES-MELBOURNE-Part1-2014-02-26.pdf
"""

columns = ['Year','US Military Fatalities per Million']

data = [[1952,190],
        [1956,  0],
        [1960,  0],
        [1964,  1],
        [1968,146],
        [1972,  0],
        [1976,  2],
        [1980,  0],
        [1984,  0],
        [1988,  0],
        [1992,  0],
        [1996,  0],
        [2000,  0],
        [2004,  4],
        [2008, 14],
        [2012,  5]]
        
deaths = pd.DataFrame(data=data, columns=columns)

# + id="5N9t0IHG123h" colab_type="code" colab={}
df = df.merge(deaths)

# + id="h1DdAls_123l" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 567} outputId="f2d07109-05d2-4ab0-ff14-f9afb14042f2"
features = ['Average Recent Growth in Personal Incomes', 
            'US Military Fatalities per Million']

target = 'Incumbent Party Vote Share'

X = df[features]
y = df[target]

model = sm.OLS(y, sm.add_constant(X))
results = model.fit()
print(results.summary())

# + id="U5JPDLhp123o" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 85} outputId="fac677a9-a5f6-458b-a6fc-8ee0cb08fdc4"
model = LinearRegression()
model.fit(X, y)
print('Intercept:', model.intercept_)
pd.Series(model.coef_, features)

# + id="wf8q2cEY123p" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 34} outputId="91a57527-256a-42b9-bab8-18b9068c508f"
np.sqrt(mean_squared_error(y, model.predict(X)))

# + [markdown] id="TovMOsNi123r" colab_type="text"
# # Train / Test Split

# + id="ku5UBe3Z123s" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 34} outputId="93724778-3151-4e5e-8d1a-974a3e5ba3ce"
train = df.query('Year < 2008')
test  = df.query('Year >= 2008')
X_train = train[features]
y_train = train[target]
X_test  = test[features]
y_test  = test[target]
X_train.shape, y_train.shape, X_test.shape, y_test.shape

# + id="mXWF1eVn123w" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 34} outputId="1cb61154-2692-4524-be1f-98dcdf8e8689"
model.fit(X_train, y_train)
model.predict(X_test)

# + id="44k2j2sv123z" colab_type="code" colab={"base_uri": "https://localhost:8080/", "height": 145} outputId="74b6faa0-e140-4626-d61e-e5a279c47988"
test

# + [markdown] id="b-8E5BX81232" colab_type="text"
# ### More about the "Bread & Peace" model
# - https://fivethirtyeight.com/features/what-do-economic-models-really-tell-us-about-elections/
# - https://statmodeling.stat.columbia.edu/2007/12/15/bread_and_peace/
# - https://avehtari.github.io/RAOS-Examples/ElectionsEconomy/hibbs.html
# - https://douglas-hibbs.com/
# - http://www.douglas-hibbs.com/HibbsArticles/HIBBS-PRESVOTE-SLIDES-MELBOURNE-Part1-2014-02-26.pdf

# + [markdown] id="aDWksVbf1233" colab_type="text"
# # Dimensionality in Linear Regression!
#
# Muliple Regression is simply an extension of the bivariate case. The reason why we see the bivariate case demonstrated so often is simply because it's easier to graph and all of the intuition from the bivariate case is the same as we keep on adding explanatory variables.
#
# As we increase the number of $x$ values in our model we are simply fitting a n-1-dimensional plane to an n-dimensional cloud of points within an n-dimensional hypercube. 

# + [markdown] id="eSErkf1J1234" colab_type="text"
# # Interpreting Coefficients
#
# One of Linear Regression's strengths is that the parameters of the model (coefficients) are readily interpretable and useful. Not only do they describe the relationship between x and y but they put a number on just how much x is associated with y. We should be careful to not speak about this relationshiop in terms of causality because these coefficients are in fact correlative measures. We would need a host of additional techniques in order to estimate a causal effect using linear regression (econometrics).
#
# \begin{align}
# \hat{\beta} = \frac{Cov(x,y)}{Var(y)}
# \end{align}
#
# Going back to the two equations for the two models that we have estimated so far, lets replace their beta values with their actual values to see if we can make sense of how to interpret these beta coefficients.

# + [markdown] id="Mrcl_weP1235" colab_type="text"
# ## Bivariate Model
#
# $y_i = \beta_0 + \beta_1temperature + \epsilon$
#
# $sales_i = -596.2 + 24.69temperature + \epsilon$
#
# What might $\beta_0$ in this model represent? It represents the level of sales that we would have if temperature were 0. Since this is negative one way of interpreting it is that it's so cold outside that you would have to pay people to eat ice cream. A more appropriate interpretation is probably that the ice cream store owner should close his store down long before the temperature reaches 0 degrees farenheit (-17.7 celsius). The owner can compare his predicted sales with his costs of doing business to know how warm the weather has to get before he should open his store.
#
# What might the $beta_1$ in this model reprsent? it represents the increase in sales for each degree of temperature increase. For every degree that the temperature goes up outside he has $25 more in sales.

# + [markdown] id="1BCj64eR1236" colab_type="text"
# ## Multiple Regression Model
#
# $y_i = \beta_0 + \beta_1age_i + \beta_2weight_i + \epsilon$
#
# $BloodPressure_i = 30.99+ .86age_i + .33weight_i + \epsilon$
#
# The interpretation of coefficients in this example are similar. The intercept value repesents the blood pressure a person would have if they were 0 years old and weighed 0 pounds. This not a super useful interpretation. If we look at our data it is unlikely that we have any measurements like these in the dataset. This means that our interpretation of our intercept likely comes from extrapolating the regression line (plane). Coefficients having straightforward interpretations is a strength of linear regression if we're careful about extrapolation and only interpreting our data within the context that it was gathered.
#
# The interpretation of our other coefficients can be a useful indicator for how much a person similar to those in our dataset's blood pressure will go up on average with each additional year of age and pound of weight.

# + [markdown] id="3X7NGUij1237" colab_type="text"
# # Basic Model Validation
#
# One of the downsides of relying on $R^2$ too much is that although it tells you when you're fitting the data well, it doesn't tell you when you're *overfitting* the data. The best way to tell if you're overfitting the data is to get some data that your model hasn't seen yet, and evaluate how your predictions do. This is essentially what "model validation" is.

# + [markdown] id="NVrAKh3S1238" colab_type="text"
# # Why is Linear Regression so Important?
#
# ## Popularity 
#
# Linear Regression is an extremely popular technique that every data scientist **needs** to understand. It's not the most advanced technique and there are supervised learning techniques that will obtain a higher accuracy, but where it lacks in accuracy it makes up for it in interpretability and simplicity.
#
# ## Interpretability
#
# Few other models possess coefficients that are so directly linked to their variables with a such a clear interpretation. Tomorrow we're going to learn about ways to make them even easier to interpret.
#
# ## Simplicity
#
# A linear regression model can be communicated just by writing out its equation. It's kind of incredible that such high dimensional relationships can be described from just a linear combination of variables and coefficients. 

StarcoderdataJupyterScriptsDedupFiltered