What function we use to estimate a relationship between \(y\) and the regressor(s) depends on the data that we have. As practice using multiple x-variables, let’s simulate a dataset that is generated by the following equation $\( y = 4 + 0.5 x + 3 x^2 + u \)$

The relationship between \(y\) and \(x\) in the above equation is said to be nonlinear in \(x\).

import pandas as pd
import numpy as np
import seaborn as sns
import statsmodels.formula.api as smf


df = pd.DataFrame(columns=['y', 'x', 'x2'])
df['x'] = np.random.normal(0, 1, 100)
df['x2'] = df['x'] ** 2
df['y'] = 4 + 0.5 * df['x'] + 3 * df['x2'] + np.random.normal(0, 0.5, 100)
y x x2
0 15.159244 1.764052 3.111881
1 4.006576 0.400157 0.160126
2 6.727911 0.978738 0.957928
3 20.669952 2.240893 5.021602
4 14.810536 1.867558 3.487773

Note that to square column 'x' to produce 'x2', all we need to do is type df2['x'] ** 2 since ** is the command for exponentiation.

Next, plot a relationship between just y and x, assuming a linear fit.

sns.lmplot(x='x', y='y', data=df)
<seaborn.axisgrid.FacetGrid at 0x7fe0343631f0>

The estimated straight line gets things horribly wrong! The line is pretty far away from a lot of the actual data points.

The scatter plot of the data above makes clear that the data demonstrates a non-linear relationship.

What we ought do to is fit a line given by $\( \hat{y} = c + \hat{\beta}_1 x + \hat{\beta}_2 x^2 \)$

where there are two beta coefficients. One for \(x\), and one for \(x^2\).

model = smf.ols('y ~ x + x2', data=df).fit()
                            OLS Regression Results                            
Dep. Variable:                      y   R-squared:                       0.983
Model:                            OLS   Adj. R-squared:                  0.983
Method:                 Least Squares   F-statistic:                     2885.
Date:                Mon, 06 Dec 2021   Prob (F-statistic):           3.91e-87
Time:                        12:06:14   Log-Likelihood:                -75.005
No. Observations:                 100   AIC:                             156.0
Df Residuals:                      97   BIC:                             163.8
Df Model:                           2                                         
Covariance Type:            nonrobust                                         
                 coef    std err          t      P>|t|      [0.025      0.975]
Intercept      4.0714      0.066     61.584      0.000       3.940       4.203
x              0.5615      0.052     10.829      0.000       0.459       0.664
x2             2.9666      0.040     73.780      0.000       2.887       3.046
Omnibus:                        7.876   Durbin-Watson:                   2.003
Prob(Omnibus):                  0.019   Jarque-Bera (JB):                3.580
Skew:                           0.188   Prob(JB):                        0.167
Kurtosis:                       2.152   Cond. No.                         2.47

[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

The variable model2 stores a lot of data. Beyond holding summary output for the performance of the regression (which we’ve accessed via the model2.summary() command), we can reference the estimated parameters directly via .params. For example:

print('Intercept: ', model.params['Intercept'])
print('beta for x:', model.params['x'])
print('beta for x2:', model.params['x2'])
Intercept:  4.071350754098632
beta for x: 0.5615154693913407
beta for x2: 2.9666242239334637

As an aside, the fact that we’re using square brackets to reference items inside of squared.params is a clue that the .params component of the variable squared was built using a dictionary-like structure.

One way to tell that the estimates of squared are better than the estimates of straight is to look at the R-squared value in the summary output. This measure takes a value between 0 and 1, with a score of 1 indicating perfect fit. For comparison purposes, consider the linear fit model below:

model_linear = smf.ols('y ~ x', data=df).fit()
                            OLS Regression Results                            
Dep. Variable:                      y   R-squared:                       0.056
Model:                            OLS   Adj. R-squared:                  0.046
Method:                 Least Squares   F-statistic:                     5.767
Date:                Mon, 06 Dec 2021   Prob (F-statistic):             0.0182
Time:                        12:06:14   Log-Likelihood:                -277.26
No. Observations:                 100   AIC:                             558.5
Df Residuals:                      98   BIC:                             563.7
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
                 coef    std err          t      P>|t|      [0.025      0.975]
Intercept      7.0734      0.392     18.054      0.000       6.296       7.851
x              0.9318      0.388      2.401      0.018       0.162       1.702
Omnibus:                       49.717   Durbin-Watson:                   1.848
Prob(Omnibus):                  0.000   Jarque-Bera (JB):              128.463
Skew:                           1.877   Prob(JB):                     1.27e-28
Kurtosis:                       7.091   Cond. No.                         1.06

[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

The linear-fit model produced an R-squared of 0.056 while the quadratic-fit model yielded a R-squared of 0.983.

Another informative way to jude a modeled relationship is by plotting the residuals. The residuals are the “un-expected” part of the equation. For example, in the linear-fit model, the residual is defined as $\( \hat{u} := y - \hat{y} = y - \hat{c} - \hat{\beta} x \)\( and in the quadratic-fit model the residual, \)\hat{u}\(, is given by \)\( \hat{u} := y - \hat{y} = y- \hat{c} - \hat{\beta}_1 x - \hat{\beta}_2 x^2 \)$ Using the information in .params, we can calculate residuals.

df['resid'] = model.resid
df['resid_linear'] = model_linear.resid

Next, plot the two sets of residuals.

sns.scatterplot(x='x', y='resid_linear', data=df)
<AxesSubplot:xlabel='x', ylabel='resid_linear'>

In the straight-line model, we can see that the errors have a noticeable pattern to them. This is an indication that a more complicated function of \(x\) would be a better description for the relationship between \(x\) and \(y\).

sns.scatterplot(x='x', y='resid', data=df)
<AxesSubplot:xlabel='x', ylabel='resid'>

In comparison, the residuals from the squared model look more random. Additionally, they’re substantially smaller on average, with almost all residuals having an absolute value less than one. This indicates a much better fit than the straight-line model in which the residual values were often much larger.

Looking Beyond OLS

OLS works well when the \(y\) variable in our model is a linear combination of \(x\) variables. Note that the relationship between \(y\) and a given regressor may be nonlinear, as in the case of \(y\) being a function of \(x\) and \(x^2\). However, while we may say that \(y\) is a nonlinear function of \(x\) in this case, the variable \(y\) is still a linear function of \(x\) and \(x^2\). To clarify: $\( y = \alpha + \beta x + u \)\( is linear in \)x\(. Likewise: \)\( y = \alpha + \beta_1 x + \beta_2 x^2 + u \)\( is linear in \)x\( and \)x^2\(. In contrast, the function \)\( y = \frac{e^{\alpha + \beta x + u}}{1 + e^{\alpha + \beta x + u}} \)\( is *nonlinear*. This last equation may look terrifyingly unnatural, but it's actually very useful. Let's get a sense of what the equation looks like by plotting the function \)\( y = \frac{e^{x}}{1 + e^{x}}. \)$

curve = pd.DataFrame(columns=['x','y'])
curve['x'] = np.random.randint(low=-10, high=10, size=100)
curve['y'] = np.exp(curve['x']) / (1 + np.exp(curve['x']))
sns.lineplot(x='x', y='y', data=curve)
    x         y
0   2  0.880797
1   5  0.993307
2 -10  0.000045
3  -7  0.000911
4  -7  0.000911
<AxesSubplot:xlabel='x', ylabel='y'>

The function \(y = e^{x}/(1+e^{x})\) creates an S-curve with a lower bound of \(0\) and an upper bound of \(1\).

These bounds are useful in analytics. Often, our task is to estimate probabilities. For instance, what is the likelihood that a borrower defaults on their mortgage, the likelihood that a credit card transaction is fraudulent, or the likelihood that a company will violate its capital expenditure covenant on an outstanding loan? All of these questions require us to estimate a probability. The above function is useful because it considers a situation in which the \(y\) variable is necessarily between \(0\) and \(1\) (i.e. \(0\%\) probability and \(100\%\) probability).

Suppose that we instead took the curve data above and tried to fit it with linear regression.

sns.lmplot(x='x', y='y', data=curve)
<seaborn.axisgrid.FacetGrid at 0x7fe005c01ac0>

Note that in the above plot, there are estimated “probabilities” (the \(\hat{y}\) values) that are either lower than \(0\) or higher than \(1\). This is statistically impossible.

That fancy S-shaped function above is called the inverse logit function. That is because \(e^x / (1+e^x)\) is the inverse of the logit function given by \(log(x / (1-x)\). Just like we can invert the equation \(y = f(x)\) to get \(f^{-1}(y) = x\), the inverse of \(y = e^x / (1+e^x)\) is given by \(log(y / (1-y)) = x\).

When we model probabilities, the \(y\) variable will be either \(0\) or \(1\) for any observations. Observations where the event occurred are recorded with \(y=1\). For instance, in a dataset about mortgage default, those borrowers who default on their mortgage would have \(y=1\) and everyone else would have \(y=0\).

Suppose that we estimate a model given by $\( log\Big(\frac{y}{1-y}\Big) = \alpha + \beta x + u \)\(. Then, the estimated value \)\hat{\alpha} + \hat{\beta}x\( for a given observation would be \)log(\hat{y}/(1-\hat{y}))\(. This is what we refer to as the *log odds ratio*. The odds ratio is the probability that \)y\( equals \)1\( divided by the probability that \)y\( equals \)0\(; this is what \)\hat{y}/(1-\hat{y})\( tells us. Ultimately, our estimate for \)\hat{y}\( then tells us the likelihood that the true value for \)y\( is equal to \)1$.

This type of model is called a logistic regression model. Let’s simulate some data.

df2 = pd.DataFrame(columns=['x', 'y'])
df2['x'] = np.random.normal(2, 1, 1000)
df2['xb'] = -9 + 4 * df2['x']
df2['y'] = np.random.binomial(n=1, p= np.exp(df2['xb']) / (1+np.exp(df2['xb'])) )
x y xb
0 3.764052 1 6.056209
1 2.400157 0 0.600629
2 2.978738 1 2.914952
3 4.240893 1 7.963573
4 3.867558 1 6.470232

Now try fitting a model using linear regression (ordinary least squares). When OLS is applied to a \(y\) variable that only takes value 0 or 1, the model is referred to as a linear probability model (LPM).

lpm = smf.ols('y ~ x', data=df2).fit()
                            OLS Regression Results                            
Dep. Variable:                      y   R-squared:                       0.514
Model:                            OLS   Adj. R-squared:                  0.514
Method:                 Least Squares   F-statistic:                     1056.
Date:                Mon, 06 Dec 2021   Prob (F-statistic):          1.41e-158
Time:                        12:06:16   Log-Likelihood:                -346.17
No. Observations:                1000   AIC:                             696.3
Df Residuals:                     998   BIC:                             706.2
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
                 coef    std err          t      P>|t|      [0.025      0.975]
Intercept     -0.2928      0.024    -12.187      0.000      -0.340      -0.246
x              0.3564      0.011     32.493      0.000       0.335       0.378
Omnibus:                      113.684   Durbin-Watson:                   2.045
Prob(Omnibus):                  0.000   Jarque-Bera (JB):               37.844
Skew:                           0.207   Prob(JB):                     6.06e-09
Kurtosis:                       2.142   Cond. No.                         5.70

[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

Now fit the model using logistic regression.

logit = smf.logit('y ~ x', data=df2).fit()
Optimization terminated successfully.
         Current function value: 0.292801
         Iterations 8
                           Logit Regression Results                           
Dep. Variable:                      y   No. Observations:                 1000
Model:                          Logit   Df Residuals:                      998
Method:                           MLE   Df Model:                            1
Date:                Mon, 06 Dec 2021   Pseudo R-squ.:                  0.5660
Time:                        12:06:16   Log-Likelihood:                -292.80
converged:                       True   LL-Null:                       -674.60
Covariance Type:            nonrobust   LLR p-value:                4.432e-168
                 coef    std err          z      P>|z|      [0.025      0.975]
Intercept     -8.5975      0.567    -15.154      0.000      -9.709      -7.486
x              3.8987      0.259     15.062      0.000       3.391       4.406

The estimated parameters (intercept and \(\beta\) coefficient) are much closer to their true value using logistic regression. While popular, linear probability models do not often give meaningful estimates. This is especially true when outcomes are rare (meaning most \(y\) values are either 0 or 1, and there is not a relatively even balance between the two).