<jupyter_start><jupyter_text>Multi-linear regression: how many variables?[![Latest release](https://badgen.net/github/release/Naereen/Strapdown.js)](https://github.com/eabarnes1010/course_objective_analysis/tree/main/code)[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/eabarnes1010/course_objective_analysis/blob/main/code/minimum_corr_for_added_value.ipynb)If I have two predictors $x_1$ and $x_2$, under what circumstances is the second one useful for predicting $y$?<jupyter_code>#.............................................
# IMPORT STATEMENTS
#.............................................
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import importlib
from sklearn import linear_model
from sklearn import metrics

mpl.rcParams['figure.facecolor'] = 'white'
mpl.rcParams['figure.dpi']= 150
dpiFig = 300.


np.random.seed(300)<jupyter_output><empty_output><jupyter_text>Let's start by creating two predictors, x1 and x2, and predictand y. x1 will be totally random, and the others will build upon that.<jupyter_code>x1 = np.random.normal(0.,1.,size=100,)
print(np.shape(x1))<jupyter_output>(100,)<jupyter_text>Now we create x2.<jupyter_code>a = 0.8
b = np.sqrt(1. - a**2)
x2 = []

# create red-noise time series iteratively
for it in np.arange(0,100,1):
    x2.append(a*x1[it] + b*np.random.normal(size=1))

x2 = np.asarray(x2)[:,0]    
print(np.shape(x2))<jupyter_output>(100,)<jupyter_text>Now let's make $y$, which is composed of pieces of x1, x2 and noise.<jupyter_code>a = 0.3
b = np.sqrt(1. - a**2)
y = []

# create red-noise time series iteratively
for it in np.arange(0,100,1):
    y.append(a*x1[it] + (.05)*x2[it] + b*np.random.normal(size=1))

y = np.asarray(y)[:,0]    
print(np.shape(y))<jupyter_output>(100,)<jupyter_text>We can calculate the correlations of the predictors and predictands just to confirm that they all have some relationship with one another.<jupyter_code>c12 = np.corrcoef(x1,x2)[0,1]
c1y = np.corrcoef(x1,y)[0,1]
c2y = np.corrcoef(y,x2)[0,1]

print('corr(x1,x2) = ' + str(np.round(c12,3)))
print('corr(x1,y)  = ' + str(np.round(c1y,3)))
print('corr(x2,y)  = ' + str(np.round(c2y,3)))<jupyter_output>corr(x1,x2) = 0.828
corr(x1,y)  = 0.426
corr(x2,y)  = 0.35<jupyter_text>Theory Based on theory, the minimum useful correlation of c2y is the following (from theory)...<jupyter_code>minUseful = np.abs(c1y*c12)
print('minimum useful corr(x2,y) = ' + str(np.round(minUseful,3)))<jupyter_output>minimum useful corr(x2,y) = 0.353<jupyter_text>Furthermore, we can show analytically that the variance explained between using x1 versus x1 and x2 is practically identical since x2 doesn't appear to add additional information (i.e. |c2y| < minUseful).<jupyter_code>#just using x1
R2 = c1y**2
print('theory: y variance explained by x1        = ' + str(np.round(R2,3)))

#using x1 and x2
R2 = (c1y**2 + c2y**2 - 2*c1y*c2y*c12)/(1-c12**2)
print('theory: y variance explained by x1 & x2   = ' + str(np.round(R2,3)))<jupyter_output>theory: y variance explained by x1        = 0.181
theory: y variance explained by x1 & x2   = 0.181<jupyter_text>Actual fits We can confirm the theory now through some fun examples where we actually fit y using x1 and x2. In fact, we see that the fits indeed give us exactly what is expected by theory.<jupyter_code># only x1 predictor
X = np.swapaxes([x1],1,0)
Y = np.swapaxes([y],1,0)

# with sklearn
regr = linear_model.LinearRegression()
regr.fit(X, Y)
R2_x1 = metrics.r2_score(Y,regr.predict(X))
print('y variance explained by x1 fit      = ' + str(np.round(R2_x1,5)))

#---------------------------------------------
# both x1 and x2 predictors
X = np.swapaxes([x1,x2],1,0)
Y = np.swapaxes([y],1,0)

# with sklearn
regr = linear_model.LinearRegression()
regr.fit(X, Y)
R2_x12 = metrics.r2_score(Y,regr.predict(X))
print('y variance explained by x1 & x2 fit = ' + str(np.round(R2_x12,5)))<jupyter_output>y variance explained by x1 fit      = 0.18142
y variance explained by x1 & x2 fit = 0.18144<jupyter_text>But what is going on here? Why is the $R^2$ slightly higher when we added x2? I thought theory said it shouldn't improve my variance explained _at all_?  What about more predictors? (aka _overfitting_)<jupyter_code>X = np.random.normal(0.,1.,size=(100,40))
Y = np.random.normal(0.,1.,size=100,)

rval = []
for n in np.arange(0,np.shape(X)[1]):
    # with sklearn
    regr = linear_model.LinearRegression()
    regr.fit(X[:,0:n+1], Y)
    R2 = metrics.r2_score(Y,regr.predict(X[:,0:n+1]))
    rval.append(R2)
    
plt.figure(figsize=(8,6))
plt.plot(np.arange(0,np.shape(X)[1]),rval,'o-')
plt.xlabel('number of random predictors')
plt.ylabel('fraction variance explained')
plt.title('Variance Explained')
plt.show()<jupyter_output><empty_output><jupyter_text>Adjusted R$^2$ There is a great solution to this - known as the _adjusted $R^2$_. It is a measure of explained variance, but you are penalized (the number decreases) when too many predictors are used. The adjusted $R^2$ increases only if the new term improves the model more than would be expected by chance.<jupyter_code>def adjustRsquared(r2,n,p):
    adjustR2 = 1 - (1-r2)*(n-1)/(n-p-1)
    return adjustR2
# only fitting with x1
p=1
n = len(x1)
adjustR2 = adjustRsquared(R2_x1,n,p)
print('fit with x1 only')
print('    R-squared = ' + str(np.round(R2_x1,3)) + ', Adjusted R-squared = ' + str(np.round(adjustR2,3)))

# fitting with x1 and x2
p = 2
n = len(x1)
adjustR2 = adjustRsquared(R2_x12,n,p)
print('fit with x1 and x2 only')
print('    R-squared = ' + str(np.round(R2_x12,3)) + ', Adjusted R-squared = ' + str(np.round(adjustR2,3)))<jupyter_output>fit with x1 only
    R-squared = 0.181, Adjusted R-squared = 0.173
fit with x1 and x2 only
    R-squared = 0.181, Adjusted R-squared = 0.165<jupyter_text>In our silly example above with 40 predictors, the adjusted R2 is the following...<jupyter_code>n = len(Y)
p = np.arange(0,np.shape(X)[1]) + 1
adjustR2 = adjustRsquared(np.asarray(rval),n,p)

plt.figure(figsize=(8,6))
plt.axhline(y=0,color='gray')
plt.plot(np.arange(1,np.shape(X)[1]+1),rval,'o-', label='R2')
plt.plot(np.arange(1,np.shape(X)[1]+1),adjustR2,'o-',color='red', label='adjusted R2')
plt.xlabel('number of predictors')
plt.ylabel('fraction variance explained')
plt.legend()
plt.title('Adjusted R-squared')
plt.show()<jupyter_output><empty_output><jupyter_text>Significance of Adjusted $R^2$To end, let's compute the adjusted R-squared many times for a lot of random data to get a feeling of the spread of possible adjusted R-squared values by chance alone.<jupyter_code>rVec = np.zeros(shape=(40,500))
for nvar in (np.arange(1,np.shape(rVec)[0]+1)):
    r = []
    for n in np.arange(0,500):
        X = np.random.normal(0.,1.,size=(100,nvar))
        Y = np.random.normal(0.,1.,size=100,)

        # with sklearn
        regr = linear_model.LinearRegression()
        regr.fit(X[:,0:n+1], Y)
        R2 = metrics.r2_score(Y,regr.predict(X[:,0:n+1]))
        r.append(R2)

    rVec[nvar-1,:] = adjustRsquared(np.asarray(r),100,nvar)
pTop = np.percentile(rVec,97.5,axis=1)
pBot = np.percentile(rVec,2.5,axis=1)

plt.figure(figsize=(8,6))
plt.axhline(y=0,color='gray')
plt.plot(np.arange(1,np.shape(X)[1]+1),adjustR2,'o-',color='red', label='adjusted R2')
plt.fill_between(np.arange(1,len(p)+1), pBot, pTop,color='lightgray', label='confidence bounds')

plt.xlabel('number of predictors')
plt.ylabel('fraction variance explained')
plt.legend()
plt.title('Adjusted R2')
plt.ylim(-1,1)
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1728509177340552, 'fraction_numerical': 0.0492539284299485, 'mean_word_length': 4.62286562731997, 'pattern_counts': {'":': 0, '<': 42, '<?xml version=': 0, '>': 41, 'https://': 4, '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': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '4551', 'license': 'MIT', 'n_tokens_mistral': 3155, 'n_tokens_neox': 2918, 'n_words': 662, 'path': 'code/minimum_corr_for_added_value.ipynb', 'repo_name': 'eabarnes1010/course_objective_analysis', 'size': '7573'}