<jupyter_start><jupyter_text>Linear regression gets to the heart ofstatistics: Given a set of data points,what is the relationship of the data in-hand to data yet seen? How should information from one data set propagate toother data? Linearregression offers the following model to address thisquestion:$$\mathbb{E}(Y|X=x) \approx a x + b$$ That is, given specific values for $X$, assume that the conditionalexpectationis a linear function of those specific values.  However, becausethe observedvalues are not the expectations themselves, the model accommodatesthis with anadditive noise term. In other words, the observed variable (a.k.a.response,target, dependent variable) is modeled as,$$\mathbb{E}(Y| X=x_i) + \epsilon_i \approx a x + b+ \epsilon_i=y$$ where $\mathbb{E}(\epsilon_i)=0$ and the $\epsilon_i$ are iid andwhere thedistribution function of $\epsilon_i$ depends on the problem, eventhough it isoften assumed Gaussian. The $X=x$ values are known as independentvariables,covariates, or regressors. Let's see if we can use all of the methods we havedeveloped so far tounderstand this form of regression. The first task is todetermine how toestimate the unknown linear parameters, $a$ and $b$. To makethis concrete,let's assume that $\epsilon \sim \mathcal{N}(0,\sigma^2)$.  Bearin mind that$\mathbb{E}(Y|X=x)$ is a deterministic function of $x$.  In otherwords, thevariable $x$ changes with each draw, but after the data have beencollectedthese are no longer random quantities. Thus, for fixed $x$, $y$ is arandomvariable generated by $\epsilon$.  Perhaps we should denote $\epsilon$ as$\epsilon_x$ to emphasize this, but because $\epsilon$ is an independent,identically-distributed (iid) random variable at each fixed $x$, this would beexcessive.  Because of Gaussian additive noise, thedistribution of $y$ iscompletely characterized by its mean and variance.$$\begin{align*}\mathbb{E}(y) &= a x  + b  \\\\mathbb{V}(y) &= \sigma^2\end{align*}$$ Using the maximum likelihood procedure, we writeout the log-likelihoodfunction as$$\mathcal{L}(a,b)  = \sum_{i=1}^n \log \mathcal{N}(a x_i +b , \sigma^2)\propto \frac{1}{2 \sigma^2}\sum_{i=1}^n (y_i-a x_i-b)^2$$ Note that we suppressed  the terms that are irrelevent tothe maximum-finding.Taking the derivative of this with respect to $a$ givesthe following equation:$$\frac{\partial \mathcal{L}(a,b)}{\partial a}= 2 \sum_{i=1}^n x_i (b+ a x_i-y_i) =0$$ Likewise, we do the same for the $b$ parameter$$\frac{\partial \mathcal{L}(a,b)}{\partial b}=2\sum_{i=1}^n (b+a x_i-y_i) =0$$ The following code simulates some data and uses Numpy tools tocompute theparameters as shown,<jupyter_code>import numpy as np
a = 6;b = 1 # parameters to estimate
x = np.linspace(0,1,100)
y = a*x + np.random.randn(len(x))+b
p,var_=np.polyfit(x,y,1,cov=True) # fit data to line
y_ = np.polyval(p,x) # estimated by linear regression
%matplotlib inline

from matplotlib.pylab import subplots
fig,axs=subplots(1,2)
fig.set_size_inches((9,3.5))
_=ax =axs[0]
_=ax.plot(x,y,'o',alpha=.3,color='gray',ms=10)
_=ax.plot(x,y_,color='k',lw=3)
_=ax.set_xlabel("$x$",fontsize=22)
_=ax.set_ylabel("$y$",fontsize=22)
_=ax.set_title("Linear regression; a=%3.3g b=%3.3g"%(p[0],p[1]))
_=ax = axs[1]
_=ax.hist(y_-y,color='gray')
_=ax.set_xlabel(r"$\Delta y$",fontsize=22)
fig.tight_layout()
fig.savefig('fig-statistics/Regression_001.png')<jupyter_output><empty_output><jupyter_text><!--dom:FIGURE: [fig-statistics/Regression_001.png, width=500 frac=1.] The panel onthe left shows the data and regression line. The panel on the right shows ahistogram of the regression errors.  --><!--begin figure -->The panel on the leftshows the data and regression line. The panel on the right shows a histogram ofthe regression errors.<img src="fig-statistics/Regression_001.png"width=500> The graph on the left of[Figure](fig:Regression_001)shows the regression line plotted against thedata. The estimatedparameters are noted in the title.  The histogram on theright of[Figure](fig:Regression_001) shows the residual errors in the model.It is always a good idea to inspect the residuals of any regressionfornormality. These are the differences between the fitted line foreach $x_i$value and the corresponding $y_i$ value in the data.Note that the $x$ term doesnot have to be uniformly monotone.To decouple the deterministic variation fromthe random variation, we can fixthe index and write separate problems of theform$$y_i = a x_i +b + \epsilon_i$$ where $\epsilon_i \sim \mathcal{N}(0,\sigma^2)$. What could we do withjustthis one component of the problem? In other words, suppose we had$m$-samples ofthis component as in $\lbrace y_{i,k}\rbrace_{k=1}^m$. Followingthe usualprocedure, we could obtain estimates of the mean of $y_i$ as$$\hat{y_i} = \frac{1}{m}\sum_{k=1}^m y_{i,k}$$ However, this tells us nothing about the individual parameters $a$and $b$because they are not separable in the terms that are computed, namely,we mayhave$$\mathbb{E}(y_i) = a x_i +b$$ but we still only have one equation and the two unknowns, $a$ and$b$. Howabout if we consider and fix another component $j$ as in$$y_j = a x_j +b + \epsilon_i$$ Then, we have$$\mathbb{E}(y_j) = a x_j +b$$ so at least now we have two equations and two unknowns and we know how toestimate the left hand sides of these equations from the data using theestimators $\hat{y_i}$ and  $\hat{y_j}$. Let's see how this works in the codesample below.<jupyter_code>x0, xn =x[0],x[80]
# generate synthetic data
y_0 = a*x0 + np.random.randn(20)+b
y_1 = a*xn + np.random.randn(20)+b
# mean along sample dimension
yhat = np.array([y_0,y_1]).mean(axis=1)
a_,b_=np.linalg.solve(np.array([[x0,1],
                                [xn,1]]),yhat)<jupyter_output><empty_output><jupyter_text><!-- dom:FIGURE: [fig-statistics/Regression_002.png, width=500 frac=0.95] Thefitted and true lines are plotted with the data values. The squares at eitherend of the solid line show the mean value for each of the data groups shown. --><divid="fig:Regression_002">The fitted and true lines are plotted withthe data values. The squares at either end of the solid line show the mean valuefor each of the data groups shown.<img src="fig-statistics/Regression_002.png" width=500>**ProgrammingTip.**The prior code uses the `solve` function in the Numpy `linalg` module,whichcontains the core linear algebra codes in Numpy that incorporate thebattle-tested LAPACK library. We can write out the solution for theestimated parameters for thiscase where $x_0 =0$$$\begin{align*}\hat{a} &= \frac{\hat{y}_i - \hat{y}_0}{x_i} \\\\hat{b} &=\hat{y_0} \\\\end{align*}$$ The expectations and variances of these estimators are the following,$$\begin{align*}\mathbb{E}(\hat{a}) &= \frac{a x_i }{x_i}=a  \\\\mathbb{E}(\hat{b}) &=b \\\\mathbb{V}(\hat{a}) &= \frac{2 \sigma^2}{x_i^2} \\\\mathbb{V}(\hat{b}) &= \sigma^2\end{align*}$$ The expectations show that the estimators are unbiased.  Theestimator$\hat{a}$ has a variance that decreases as larger points $x_i$ areselected.That is, it is better to have samples further out along thehorizontal axis forfitting the line.  This variance quantifies the *leverage*of those distantpoints.**Regression From Projection Methods.** Let's see if we can apply ourknowledgeof projection methods to the general case.  In vector notation, we canwritethe following:$$\mathbf{y} = a \mathbf{x} + b\mathbf{1} + \boldsymbol{\epsilon}$$ where $\mathbf{1}$ is the vector of all ones. Let's use the inner-productnotation,$$\langle \mathbf{x},\mathbf{y} \rangle = \mathbb{E}(\mathbf{x}^T \mathbf{y})$$ Then, by taking the inner-product with some $\mathbf{x}_1 \in\mathbf{1}^\perp$we obtain [^perp],[^perp]: The space of all vectors, $\mathbf{a}$ such that$\langle\mathbf{a},\mathbf{1} \rangle = 0$ is denoted $\mathbf{1}^\perp$.$$\langle \mathbf{y},\mathbf{x}_1  \rangle = a \langle \mathbf{x},\mathbf{x}_1\rangle$$ Recall that $\mathbb{E}(\boldsymbol{\epsilon})=\mathbf{0}$.  Wecan finallysolve for $a$ as$$\begin{equation}\hat{a} = \frac{\langle\mathbf{y},\mathbf{x}_1 \rangle}{\langle\mathbf{x},\mathbf{x}_1 \rangle} \end{equation}\label{eq:ahat} \tag{1}$$ That was pretty neat but now we have the mysterious $\mathbf{x}_1$vector.Where does this come from?  If we project $\mathbf{x}$ onto the$\mathbf{1}^\perp$, then we get the MMSE  approximation to $\mathbf{x}$ in the$\mathbf{1}^\perp$ space. Thus, we take$$\mathbf{x}_1 = P_{\mathbf{1}^\perp} (\mathbf{x})$$ Remember that $P_{\mathbf{1}^\perp} $ is a projection matrix so the length of$\mathbf{x}_1$ is at most $\mathbf{x}$. This means that the denominator inthe$\hat{a}$ equation above is really just the length of the $\mathbf{x}$vector inthe coordinate system of $P_{\mathbf{1}^\perp} $. Because theprojection isorthogonal (namely, of minimum length), the Pythagorean theoremgives thislength as the following:$$\langle \mathbf{x},\mathbf{x}_1 \rangle ^2=\langle \mathbf{x},\mathbf{x}\rangle- \langle\mathbf{1},\mathbf{x} \rangle^2$$ The first term on the right is the length of the $\mathbf{x}$ vectorand lastterm is the length of $\mathbf{x}$ in the coordinate system orthogonalto$P_{\mathbf{1}^\perp} $, namely that of $\mathbf{1}$.  Wecan use this geometricinterpretation to understand what is going on intypical linear regression inmuch more detail. The fact that the denominator isthe orthogonal projection of$\mathbf{x}$ tells us that the choice of$\mathbf{x}_1$ has the strongest effect(i.e., largest value) on reducingthe variance of $\hat{a}$. That is, the more$\mathbf{x}$ is aligned with$\mathbf{1}$, the worse the variance of $\hat{a}$.This makes intuitive sensebecause the closer $\mathbf{x}$ is to $\mathbf{1}$,the more constant it is,and we have already seen from our one-dimensionalexample that distance betweenthe $x$ terms pays off in reduced variance. Wealready know that $\hat{a}$is an unbiased estimator, and, because we chose$\mathbf{x}_1$ deliberately as aprojection, we know that it is also of minimumvariance. Such estimators areknown as Minimum-Variance Unbiased Estimators(MVUE).In the same spirit, let's examine the numerator of $\hat{a}$ inEquation [1](eq:ahat). We can write$\mathbf{x}_{1}$ as the following$$\mathbf{x}_{1} = \mathbf{x} -  P_{\mathbf{1}} \mathbf{x}$$ where $P_{\mathbf{1}}$ is projection matrix  of $\mathbf{x}$ onto the$\mathbf{1}$ vector. Using this, the numerator of $\hat{a}$ becomes$$\langle \mathbf{y}, \mathbf{x}_1\rangle  =\langle \mathbf{y},\mathbf{x}\rangle -\langle \mathbf{y}, P_{\mathbf{1}} \mathbf{x}\rangle$$ Note that,$$P_{\mathbf{1}}  = \mathbf{1} \mathbf{1}^T \frac{1}{n}$$ so that writing this out explicitly gives$$\langle \mathbf{y}, P_{\mathbf{1}} \mathbf{x}\rangle = \left(\mathbf{y}^T\mathbf{1}\right) \left(\mathbf{1}^T \mathbf{x}\right)/n =\left(\sumy_i\right)\left(\sum x_{i}\right)/n$$ and similarly, we have the following for the denominator:$$\langle \mathbf{x}, P_{\mathbf{1}} \mathbf{x}\rangle = \left(\mathbf{x}^T\mathbf{1}\right) \left(\mathbf{1}^T \mathbf{x}\right)/n =\left(\sumx_i\right)\left(\sum x_{i}\right)/n$$ So, plugging all of this together gives the following,$$\hat{a} = \frac{\mathbf{x}^T\mathbf{y}-(\sum x_i)(\sum y_i)/n}{\mathbf{x}^T\mathbf{x}  -(\sum x_i)^2/n}$$ with corresponding variance,$$\begin{align*}\mathbb{V}(\hat{a}) &= \sigma^2\frac{\|\mathbf{x}_1\|^2}{\langle\mathbf{x},\mathbf{x}_1\rangle^2} \\\&= \frac{\sigma^2}{\Vert \mathbf{x}\Vert^2-n(\overline{x}^2)} \end{align*}$$ Using the same approach with $\hat{b}$ gives,$$\begin{equation}\hat{b}  = \frac{\langle \mathbf{y},\mathbf{x}^{\perp}\rangle}{\langle \mathbf{1},\mathbf{x}^{\perp}\rangle} \label{_auto1} \tag{2}\end{equation}$$$$\begin{equation} \         = \frac{\langle\mathbf{y},\mathbf{1}-P_{\mathbf{x}}(\mathbf{1})\rangle}{\langle\mathbf{1},\mathbf{1}-P_{\mathbf{x}}(\mathbf{1})\rangle} \label{_auto2} \tag{3}\end{equation}$$$$\begin{equation} \         = \frac{\mathbf{x}^T \mathbf{x}(\sum y_i)/n-\mathbf{x}^T\mathbf{y}(\sum x_i)/n}{\mathbf{x}^T \mathbf{x} -(\sum x_i)^2/n}\label{_auto3} \tag{4}\end{equation}$$ where$$P_{\mathbf{x}} = \frac{\mathbf{\mathbf{x} \mathbf{x}^T}}{\| \mathbf{x} \|^2}$$ with variance$$\begin{align*}\mathbb{V}(\hat{b})&=\sigma^2 \frac{\langle\boldsymbol{\mathbf{1}-P_{\mathbf{x}}(\mathbf{1})},\boldsymbol{\mathbf{1}-P_{\mathbf{x}}(\mathbf{1})}\rangle}{\langle\mathbf{1},\boldsymbol{\mathbf{1}-P_{\mathbf{x}}(\mathbf{1})}\rangle^2} \\\&=\frac{\sigma^2}{n-\frac{(n\overline{x})^2}{\Vert\mathbf{x}\Vert^2}}\end{align*}$$**Qualifying the Estimates.** Our formulas for the variance above include theunknown $\sigma^2$, which we must estimate from the data itself using ourplug-in estimates.  We can form the residual sum of squares as$$\texttt{RSS} = \sum_i (\hat{a} x_i + \hat{b} - y_i)^2$$ Thus, the estimate of $\sigma^2$ can be expressed as$$\hat{\sigma}^2 = \frac{\texttt{RSS}}{n-2}$$ where $n$ is the number of samples. This is also known as the*residual meansquare*. The $n-2$  represents the *degrees of freedom* (`df`).Because weestimated two parameters from the same data we have $n-2$ instead of$n$. Thus,in general, $\texttt{df} = n - p$, where $p$ is the number ofestimatedparameters.  Under the assumption that the noise is Gaussian,  the$\texttt{RSS}/\sigma^2$ is chi-squared distributed with $n-2$ degrees offreedom. Another important term is the *sum of squares about the mean*, (a.k.a*corrected* sum of squares),$$\texttt{SYY} = \sum (y_i - \bar{y})^2$$ The $\texttt{SYY}$ captures the idea of not using the $x_i$ data andjust usingthe mean of the $y_i$ data to estimate $y$. These two terms leadto the $R^2$term,$$R^2=1-\frac{\texttt{RSS}}{ \texttt{SYY} }$$ Note that for perfect regression, $R^2=1$.  That is, if theregression getseach $y_i$ data point exactly right, then$\texttt{RSS}=0$ this term equals one.Thus, this term is used tomeasure of goodness-of-fit. The `stats` module in`scipy` computesmany of these terms automatically,<jupyter_code>from scipy import stats
slope,intercept,r_value,p_value,stderr = stats.linregress(x,y)<jupyter_output><empty_output><jupyter_text>where the square of the `r_value` variable is the $R^2$ above. Thecomputedp-value is the two-sided hypothesis test with a null hypothesis thatthe slopeof the line is zero. In other words, this tests whether or not thelinearregression  makes sense for the data for that hypothesis.  TheStatsmodelsmodule provides a powerful extension to Scipy's stats module bymaking it easyto do regression and keep track of these parameters. Let'sreformulate ourproblem using the Statsmodels framework by creatinga Pandas dataframe for thedata,<jupyter_code>import statsmodels.formula.api as smf
from pandas import DataFrame
import numpy as np
d = DataFrame({'x':np.linspace(0,1,10)}) # create data
d['y'] = a*d.x+ b + np.random.randn(*d.x.shape)<jupyter_output><empty_output><jupyter_text>Now that we have  the input data in the abovePandas dataframe, wecan perform the regression as in the following,<jupyter_code>results = smf.ols('y ~ x', data=d).fit()<jupyter_output><empty_output><jupyter_text>The $\sim$ symbol is notation for $y = a x + b + \epsilon$, where theconstant$b$ is implicit in this usage of Statsmodels. The names in the stringare takenfrom the columns in the dataframe. This makes it very easy to buildmodels withcomplicated interactions between the named columns in thedataframe. We canexamine a report of the model fit by looking at the summary,<jupyter_code>print (results.summary2())<jupyter_output>Results: Ordinary least squares
=================================================================
Model:              OLS              Adj. R-squared:     0.780   
Dependent Variable: y                AIC:                32.6081 
Date:               2019-08-27 12:16 BIC:                33.2132 
No. Observations:   10               Log-Likelihood:     -14.304 
Df Model:           1                F-statistic:        32.89   
Df Residuals:       8                Prob (F-statistic): 0.000436
R-squared:          0.804            Scale:              1.2790  
-------------------------------------------------------------------
             Coef.    Std.Err.     t      P>|t|     [0.025   0.975]
-------------------------------------------------------------------
Intercept    1.4044     0.6647   2.1128   0.0676   -0.1284   2.9372
x            6.4271     1.1206   5.7354   0.0004    3.8430   9.0112
-----------------------------------------------------------------
Omnibus:         [...]<jupyter_text>```Results: Ordinary least squares=================================================================Model:              OLS              Adj. R-squared:     0.808   Dependent Variable: y                AIC:                28.1821 Date:               0000-00-00 00:00 BIC:                00.0000 No. Observations:   10               Log-Likelihood:     -12.091 Df Model:           1                F-statistic:        38.86   Df Residuals:       8                Prob (F-statistic): 0.000250R-squared:          0.829            Scale:              0.82158 -------------------------------------------------------------------              Coef.    Std.Err.     t      P>|t|    [0.025   0.975]-------------------------------------------------------------------Intercept     1.5352     0.5327   2.8817   0.0205   0.3067   2.7637x             5.5990     0.8981   6.2340   0.0003   3.5279   7.6701There is a lot more here than we have discussed so far, but theStatsmodelsdocumentation is the best place to go for complete informationabout thisreport. The F-statistic attempts to capture the contrast betweenincluding theslope parameter or leaving it off. That is, consider twohypotheses:```$$\begin{align*}H_0 \colon  \mathbb{E}(Y|X=x) &= b \\\H_1 \colon\mathbb{E}(Y|X=x) &= b + a x\end{align*}$$ In order to quantify how much better adding the slope term is fortheregression, we compute  the following:$$F = \frac{\texttt{SYY} - \texttt{RSS}}{ \hat{\sigma}^2 }$$ The numerator computes the difference in the residual squared errorsbetweenincluding the slope in the regression or just using the mean of the$y_i$values. Once again, if we assume (or can claim asymptotically) that the$\epsilon$ noise term is Gaussian, $\epsilon \sim \mathcal{N}(0,\sigma^2)$,thenthe $H_0$ hypothesis will follow an F-distribution [^fdist] with degrees offreedom from the numerator and denominator. In this case, $F \sim F(1,n-2)$.Thevalue of this statistic is reported by Statsmodels above. The correspondingreported probability shows the chance of $F$ exceeding its computed value  if$H_0$ were true. So, the take-home message from all this is that including theslope leads to a much smaller reduction in squared error than could be expectedfrom a favorable draw of $n$ points of this data, under the Gaussian additivenoise assumption. This is evidence that including the slope is meaningful forthis data.[^fdist]: The $F(m,n)$ F-distribution has two integer degree-of-freedom parameters, $m$and $n$.The Statsmodels report also shows theadjusted $R^2$ term.This is a correction to the $R^2$ calculation that accountsfor the number of parameters $p$ that the regression isfitting and the samplesize $n$,$$\texttt{Adjusted } R^2 = 1- \frac{\texttt{RSS}/(n-p)}{\texttt{SYY}/(n-1)}$$ This is always lower than $R^2$ except when $p=1$ (i.e., estimatingonly $b$).This becomes a better way to compare regressions when one isattempting to fitmany parameters with comparatively small $n$.**Linear Prediction.**  Usinglinear regression for prediction introducessome other issues. Recall thefollowing expectation,$$\mathbb{E}(Y|X=x) \approx \hat{a} x + \hat{b}$$ where we have determined $\hat{a}$  and $\hat{b}$  from the data.Given a newpoint of interest, $x_p$, we would certainly compute$$\hat{y}_p = \hat{a} x_p + \hat{b}$$ as the predicted value for $\hat{ y_p }$. This is the same as sayingthat ourbest prediction for $y$ based on $x_p$ is the above conditionalexpectation. Thevariance for this is the following,$$\mathbb{V}(y_p) = x_p^2 \mathbb{V}(\hat{a}) +\mathbb{V}(\hat{b})+2 x_p\texttt{cov}(\hat{a}\hat{b})$$ Note that we have the covariance above because $\hat{a}$ and$\hat{b}$ arederived from the same data.  We can work this out below usingour previousnotation from [1](eq:ahat),$$\begin{align*}\texttt{cov}(\hat{a}\hat{b})=&\frac{\mathbf{x}_1^T\mathbb{V}\lbrace\mathbf{y}\mathbf{y}^T\rbrace\mathbf{x}^{\perp}}{(\mathbf{x}_1^T\mathbf{x})(\mathbf{1}^T \mathbf{x}^{\perp})} = \frac{\mathbf{x}_1^T\sigma^2\mathbf{I}\mathbf{x}^{\perp}}{(\mathbf{x}_1^T \mathbf{x})(\mathbf{1}^T\mathbf{x}^{\perp})}\\\=&\sigma^2\frac{\mathbf{x}_1^T\mathbf{x}^{\perp}}{(\mathbf{x}_1^T\mathbf{x})(\mathbf{1}^T \mathbf{x}^{\perp})} =\sigma^2\frac{\left(\mathbf{x}-P_1\mathbf{x}\right)^T\mathbf{x}^{\perp}}{(\mathbf{x}_1^T\mathbf{x})(\mathbf{1}^T \mathbf{x}^{\perp})}\\\=&\sigma^2\frac{-\mathbf{x}^T P_1^T\mathbf{x}^{\perp}}{(\mathbf{x}_1^T\mathbf{x})(\mathbf{1}^T \mathbf{x}^{\perp})} =\sigma^2\frac{-\mathbf{x}^T\frac{1}{n}\mathbf{1}\mathbf{1}^T\mathbf{x}^{\perp}}{(\mathbf{x}_1^T \mathbf{x})(\mathbf{1}^T\mathbf{x}^{\perp})}\\\=&\sigma^2\frac{-\mathbf{x}^T\frac{1}{n}\mathbf{1}}{(\mathbf{x}_1^T \mathbf{x})}= \frac{-\sigma^2\overline{x}}{\sum_{i=1}^n(x_i^2-\overline{x}^2)}\\\\end{align*}$$ After plugging all this in, we obtain the following,$$\mathbb{V}(y_p)=\sigma^2 \frac{x_p^2-2 x_p\overline{x}+\Vert\mathbf{x}\Vert^2/n}{\Vert\mathbf{x}\Vert^2-n\overline{x}^2}$$ where, in practice, we use the plug-in estimate forthe  $\sigma^2$.There isan important consequence for the confidence interval for$y_p$. We cannot simplyuse the square root of $\mathbb{V}(y_p)$to form the confidence interval becausethe model includes theextra $\epsilon$ noise term. In particular, theparameters werecomputed using a set of statistics from the data, but now mustinclude different realizations for the noise term for theprediction part.  Thismeans we have to compute$$\eta^2=\mathbb{V}(y_p)+\sigma^2$$ Then, the 95\% confidence interval $y_p \in(y_p-2\hat{\eta},y_p+2\hat{\eta})$is the following,$$\mathbb{P}(y_p-2\hat{\eta}< y_p <y_p+2\hat{\eta})\approx\mathbb{P}(-2<\mathcal{N}(0,1)<2) \approx 0.95$$ where $\hat{\eta}$ comes from substituting theplug-in estimate for $\sigma$. Extensions to Multiple CovariatesWith all the machinery we have, it is ashort notational hop to considermultiple regressors as in the following,$$\mathbf{Y} = \mathbf{X} \boldsymbol{\beta} +\boldsymbol{\epsilon}$$ with the usual $\mathbb{E}(\boldsymbol{\epsilon})=\mathbf{0}$ and$\mathbb{V}(\boldsymbol{\epsilon})=\sigma^2\mathbf{I}$. Thus, $\mathbf{X}$ is a$n \times p$ full rank matrix of regressors and $\mathbf{Y}$ is the $n$-vectorof observations. Note that the constant term has been incorporated into$\mathbf{X}$ as a column of ones.  The corresponding estimated solution for$\boldsymbol{\beta}$  is the following,$$\hat{ \boldsymbol{\beta} } = (\mathbf{X}^T \mathbf{X})^{-1} \mathbf{X}^T\mathbf{Y}$$ with corresponding variance,$$\mathbb{V}(\hat{\boldsymbol{\beta}})=\sigma^2(\mathbf{X}^T \mathbf{X})^{-1}$$ and with the assumption of Gaussian errors, we have$$\hat{\boldsymbol{\beta}}\sim \mathcal{N}(\boldsymbol{\beta},\sigma^2(\mathbf{X}^T \mathbf{X})^{-1})$$ The unbiased estimate of $\sigma^2$ is the following,$$\hat{\sigma}^2 = \frac{1}{n-p}\sum \hat{\epsilon}_i^2$$ where $\hat{ \boldsymbol{\epsilon}}=\mathbf{X}\hat{\boldsymbol{\beta}}-\mathbf{Y}$ is the vector of residuals. Tukey christened the following matrixas the *hat*matrix (a.k.a. influence matrix),$$\mathbf{V}=\mathbf{X}(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T$$ because it maps $\mathbf{Y}$ into $\hat{ \mathbf{Y} }$,$$\hat{ \mathbf{Y} } = \mathbf{V} \mathbf{Y}$$ As an exercise you can check that $\mathbf{V}$ is a projectionmatrix. Notethat that matrix is solely a function of $\mathbf{X}$. Thediagonal elements of$\mathbf{V}$ are called the *leverage values* and arecontained in the closedinterval $[1/n,1]$. These terms measure of distancebetween the values of $x_i$and the mean values over the $n$ observations.Thus, the leverage terms dependonly on $\mathbf{X}$. This is thegeneralization of our initial discussion ofleverage where we had multiplesamples at only two $x_i$ points. Using the hatmatrix, we can compute thevariance of each residual, $e_i = \hat{y}-y_i$ as$$\mathbb{V}(e_i) = \sigma^2 (1-v_{i})$$ where $v_i=V_{i,i}$. Given the above-mentioned bounds on $v_{i}$,these arealways less than $\sigma^2$.Degeneracy in the columns of $\mathbf{X}$ canbecome a problem. This is whentwo or more of the columns become co-linear. Wehave already seen this with oursingle regressor example wherein $\mathbf{x}$close to $\mathbf{1}$ was badnews. To compensate for this effect we can loadthe diagonal elements and solvefor the unknown parameters as in the following,$$\hat{ \boldsymbol{\beta} } = (\mathbf{X}^T \mathbf{X}+\alpha \mathbf{I})^{-1}\mathbf{X}^T \mathbf{Y}$$ where $\alpha>0$ is a tunable hyper-parameter. This method is knownas *ridgeregression* and was proposed in 1970 by Hoerl and Kenndard. It can beshown thatthis is the equivalent to minimizing the following objective,$$\Vert \mathbf{Y}- \mathbf{X} \boldsymbol{\beta}\Vert^2  + \alpha \Vert\boldsymbol{\beta}\Vert^2$$ In other words, the length of the estimated $\boldsymbol{\beta}$ ispenalizedwith larger $\alpha$. This has the effect of stabilizing thesubsequent inversecalculation and also providing a means to trade bias andvariance, which we willdiscuss at length in the section[ch:ml:sec:regularization](ch:ml:sec:regularization)**InterpretingResiduals.**  Our model assumes an additive Gaussian noise term.We can checkthe voracity of this assumption by examining the residuals afterfitting.  Theresiduals are the difference between the fitted values and theoriginal data$$\hat{\epsilon}_i = \hat{a} x_i + \hat{b} - y_i$$ While the p-value and the F-ratio provide some indication of whetheror notcomputing the slope of the regression makes sense, we can get directlyat thekey assumption of additive Gaussian noise.For sufficiently small dimensions,the `scipy.stats.probplot` we discussed inthe last chapter provides quickvisual evidence one way or another by plottingthe standardized residuals,$$r_i = \frac{e_i}{\hat{\sigma}\sqrt{1-v_i}}$$ The other part of the iid assumption implies homoscedasticity (all$r_i$ haveequal variances). Under the additive Gaussian noise assumption, the$e_i$ shouldalso be distributed according to $\mathcal{N}(0,\sigma^2(1-v_i))$.Thenormalized residuals $r_i$ should then be distributed according to$\mathcal{N}(0,1)$. Thus, the presence of any $r_i \notin [-1.96,1.96]$ shouldnot be common at the 5% significance level and is thereby breeds suspicionregarding the homoscedasticity assumption.The Levene test in`scipy.stats.leven` tests the null hypothesis that all thevariances are equal.This basically checks whether or not the standardizedresiduals vary across$x_i$ more than expected. Under the homoscedasticityassumption, the varianceshould be independent of $x_i$. If not, then this is aclue that there is amissing variable in the analysis or that the variablesthemselves should betransformed (e.g., using the $\log$ function) into anotherformat that canreduce this effect. Also, we can use weighted least-squaresinstead of ordinaryleast-squares.**Variable Scaling.** It is tempting to conclude in a multipleregression that small coefficients in any of the $\boldsymbol{\beta}$ termsimplies that those terms are not important.  However, simple unit conversionscan cause this effect. For example, if one of the regressors is inunits ofkilometers and the others are in meters, then just the scalefactor can give theimpression of outsized or under-sized effects. Thecommon way to account forthis is to scale the regressors so that$$x^\prime = \frac{x-\bar{x}}{\sigma_x}$$ This has the side effect of converting the slope parametersinto correlationcoefficients, which is bounded by $\pm 1$.**Influential Data.** We havealready discussed the ideaof leverage. The concept of *influence* combinesleverage withoutliers. To understand influence, consider[Figure](fig:Regression_005).<!-- dom:FIGURE: [fig-statistics/Regression_005.png, width=500 frac=0.65] The point on the right hasoutsized influence in this data because it is the only one used to determine theslope of the fitted line.   --><!-- beginfigure -->The point on the right hasoutsized influence in this data because it is the only one used to determine theslope of the fitted line.<img src="fig-statistics/Regression_005.png"width=500>The point on the right in[Figure](fig:Regression_005) is the only one thatcontributes to thecalculation of the slope for the fitted line. Thus, itis very influential inthis sense. Cook's distance is a good way to get at thisconcept numerically.To compute this, we have to compute the $j^{th}$component of the estimatedtarget variable with the $i^{th}$ point deleted. Wecall this $\hat{y}_{j(i)}$.Then, we compute the following,$$D_i =\frac{\sum_j (\hat{y}_j- \hat{y}_{j(i)})^2}{p/n \sum_j (\hat{y}_j-y_j)^2}$$ where, as before, $p$ is the number of estimated terms(e.g., $p=2$ in thebivariate case). This calculation emphasizes theeffect of the outlier bypredicting the target variable with andwithout each point. In the case of[Figure](fig:Regression_005),losing any of the points on the left cannotchange the estimatedtarget variable much, but losing the single point on theright surelydoes. The point on the right does not seem to be an outlier (it*is*on the fitted line), but this is because it is influential enough torotatethe line to align with it. Cook's distance helps capture thiseffect by leavingeach sample out and re-fitting the remainder asshown in the last equation.[Figure](fig:Regression_006) shows thecalculated Cook's distance for the datain [Figure](fig:Regression_005), showing that the data point on the right(sample index `5`) has outsized influence on the fitted line. As arule ofthumb, Cook's distance values that  are larger than one aresuspect.<!--dom:FIGURE: [fig-statistics/Regression_006.png, width=500 frac=0.65] Thecalculated Cook's distance for the data in [Figure](fig:Regression_005).  <divid="fig:Regression_006"> --><divid="fig:Regression_006">The calculated Cook's distance for the datain [Figure](fig:Regression_005).<img src="fig-statistics/Regression_006.png" width=500>As anotherillustration of influence, consider [Figure](fig:Regression_007) which showssome data that nicely line up, butwith one outlier (filled black circle) in theupper panel.  The lowerpanel shows so-computed Cook's distance for this dataand emphasizesthe presence of the outlier.  Because the calculation involvesleavinga single sample out and re-calculating the rest, it can be atime-consuming operation suitable to relatively small data sets. Thereis always thetemptation to downplay the importance of outliersbecause they conflict with afavored model, but outliers must becarefully examined to understand why themodel is unable to capturethem. It could be something as simple as faulty datacollection, or itcould be an indication of deeper issues that have beenoverlooked. Thefollowing code shows how Cook's distance was compute for[Figure](fig:Regression_006) and [Figure](fig:Regression_007).<jupyter_code>fit = lambda i,x,y: np.polyval(np.polyfit(x,y,1),i)
omit = lambda i,x: ([k for j,k in enumerate(x) if j !=i])
def cook_d(k):
    num = sum((fit(j,omit(k,x),omit(k,y))-fit(j,x,y))**2 for j in x)
    den = sum((y-np.polyval(np.polyfit(x,y,1),x))**2/len(x)*2)
    return num/den<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15169146636143965, 'fraction_numerical': 0.021111535037517486, 'mean_word_length': 4.540426281486701, 'pattern_counts': {'":': 0, '<': 45, '<?xml version=': 0, '>': 44, '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': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '353515', 'license': 'MIT', 'n_tokens_mistral': 10016, 'n_tokens_neox': 9401, 'n_words': 3360, 'path': 'chapter/statistics/Regression.ipynb', 'repo_name': 'derakding/Python-for-Probability-Statistics-and-Machine-Learning-2E', 'size': '31452'}