From Simple Linear Regression Analysis to Covariance & Correlation to Independent Determinant, and R-Squared

[Eli]

Regression analysis (see how it was recently used, among other reasons, to nullify the Kenyan general election) is a means of establishing if there is a statistically significant relationship between two sets of variables, usually a set of independent variables and a set of dependent variables. On the other hand, linear regression analysis is used to forecast new observations. This means that, we can use the relationship between the variables we already know to forecast unobserved values. The idea is to understand how the value of the dependent variable "or criterion variable", say $y$ varies as we change the value of the independent variable, say $x$ while the other values of independent variable are held fixed . The dependent variable in this case, is the variable whose values we want to explain or forecast, and its values depend on the other terms in the equation. Independent variable is the one from which we can derive values of the dependent variable.

The basis of linear regression analysis is to find a straight line (linear equation) that fits or forecasts a set of observations/unobserved values in two dimensions such as all points in a scatter plot which follow or are expected to follow some pattern (usually linear).

Suppose we want to find the relationship between student grades $y$ and the time $x$ spent in studying. The linear regression equation of the grade estimates is given by

$$\begin{array}{}\text{(1)}& \hat{y}=b_0+b_{1}x\end{array}$$
where

• $\hat{y}$ is the grade estimates,
• $b_0$ is the $y$ intercept,
• $b_1$ is the gradient/slope of the line. The slope or gradient determines the sensitivity of the dependent variable $y$ on changes of the values of independent variable $x$ - how fast/slow the dependent variable varies with an independent variable.
• $x$ is the student's study time.

Now, let's represent our observations, student's study time $x$ in hours versus grades/scores $y$ in a tabular form:


$$\begin{array}{|cc|}\hline x& y\\ 1& 2\\ 2& 4\\ 3& 5\\ 4& 4\\ 5& 5\\ \hline\end{array}$$

We use a simple Python code below to plot the regression line, together with the lines that indicate the means of independent and dependent variables.

Code: [Select all] [Expand/Collapse]

1. #Import required libraries
2. import numpy as np
3. import matplotlib.pyplot as plt
4. #Create data points
5. x = np.array([1,2,3,4,5])
6. y = np.array([2,4,5,4,5])
7.      
8. fig  = plt.subplot(111, aspect='auto')
9. plt.xlabel("Independent variable x")
10. plt.ylabel("Dependent variable  y")
11. plt.ylim(top=7)
12. fig.plot(x, y, 'o', markersize=7, color='black', alpha=1.0, label="Scatter plot")
13. plt.axvline(x=np.mean(x), color='red', linestyle='-', alpha=1.0, label="x Mean")
14. plt.axhline(y=np.mean(y), color='green', linestyle='-', alpha=1.0, label="y Mean")
15. #Fit function
16. b1, b0 = np.polyfit(x, y, deg=1)
17. f = lambda x: b1*x + b0
18. plt.plot(x,f(x), c="purple", linestyle = "-", label="Regression line")
19. #Legend, save
20. plt.legend(loc=2)
21. plt.savefig('regression_line.png', format="png", bbox_inches='tight')

It turns out that the regression line passes through the point where the lines for the means of independent variable $x$ and dependent variable $y$ cross, in our case, this point is $(3,4)$, see the figure below.



The real world observational data do not always behave linearly. In general, the data will be a series of $x$ and $y$ observations (see the table and the figure above), which on average may follow a linear pattern, making it possible to represent them by a linear regression model. However, as we have seen from the figure above, the linear regression line may not always intersect or pass through any of these $(x,y)$ data points, rather, there are going to be errors (that can be measured) due to deviation of data points from the regression line for even the largely unobserved population of values of the independent and dependent variables.

Errors and Error bars

We introduce these errors (or unobserved deviations from the regression line), and denote them as $ϵ$, which are distances between the data from the assumed real life (Scatter plot) and the actual linear regression line (purple line). We modify the code above to include and illustrate the errors.

Code: [Select all] [Expand/Collapse]

1. #Import required libraries
2. import numpy as np
3. import matplotlib.pyplot as plt
4. #Create data points
5. x = np.array([1,2,3,4,5])
6. y = np.array([2,4,5,4,5])
7.      
8. fig  = plt.subplot(111, aspect='auto')
9. plt.title(" ")
10. plt.xlabel("Independent variable x")
11. plt.ylabel("Dependent variable  y")
12. plt.ylim(top=7)
13. fig.plot(x, y, 'o', markersize=7, color='black', alpha=1.0, label="Scatter plot")
14. plt.axvline(x=np.mean(x), color='red', linestyle='--', alpha=1.0, label="x Mean")
15. plt.axhline(y=np.mean(y), color='green', linestyle='--', alpha=1.0, label="y Mean")
16. #Get intercept and the slope of the fit function
17. b1, b0 = np.polyfit(x, y, deg=1)#Fit as first degree polynomial
18. f = lambda x: b1*x + b0
19. plt.plot(x,f(x), c="purple", linestyle = "-", label="Regression line")
20. #Plot half errorbars
21. plt.errorbar(x,f(x),xerr=0, yerr=[np.zeros(5), y-f(x)], fmt="o", c="blue", label=r"${\rm Error,} \ \epsilon$")
22. plt.legend(loc=2)
23. plt.savefig('regression_line_errors.png', format="png", bbox_inches='tight')

The result is the figure below that include the errors, indicated by vertical blue lines:



You could change yerr=[np.zeros(5), y-f(x)] to yerr=[y-f(x), np.zeros(5)] to indicate the opposite errors if needed,



but a standard practice usually includes both negative and positive errors, because we only care about the magnitude of the errors, not direction. You can achieve this by replacing yerr=[np.zeros(5), y-f(x)] with yerr=y-f(x) in the code above:



Depending on your data type and the errors you're interested in, below are few illustrations of typical error bars.

Example 1: error bars in x and y

Code: [Select all] [Expand/Collapse]

1. import numpy as np
2. import matplotlib.pyplot as plt
3. #Generate xy data
4. x = np.arange(0.1, 4, 0.5)
5. y = np.exp(-x)
6.  
7. #Variable error bar values
8. yerr = 0.1 + 0.2*np.sqrt(x)
9. xerr = 0.1 + yerr
10.  
11. #Illustrate basic Matplotlib (can use Pylab) ploting interface
12. plt.figure()
13. plt.errorbar(x, y, xerr=0.2, yerr=0.4)
14. plt.title("Errorbars in x and y")
15. plt.savefig('xy_Errobars.png', format="png", bbox_inches='tight')

Code: [Select all] [Expand/Collapse]

1. import numpy as np
2. import matplotlib.pyplot as plt
3.  
4. fig = plt.figure(0)
5. x = np.arange(10.0)
6. y = np.sin(np.arange(10.0) / 20.0 * np.pi)
7.  
8. plt.errorbar(x, y,xerr=0.5, yerr=0.1)
9.  
10. y = np.sin(np.arange(10.0) / 20.0 * np.pi) + 1
11. plt.errorbar(x, y,xerr=0.4, yerr=0.1, uplims=True)
12.  
13. y = np.sin(np.arange(10.0) / 20.0 * np.pi) + 2
14.  
15. upperlimits = np.array([1, 0] * 5)
16. lowerlimits = np.array([0, 1] * 5)
17. plt.errorbar(x, y, xerr=0.5, yerr=0.1, uplims=upperlimits, lolims=lowerlimits)
18. plt.savefig('Errobars_xy.png', format="png", bbox_inches='tight')
19. plt.xlim(-1, 10)

Example 2 : x error bars only

Code: [Select all] [Expand/Collapse]

1. import numpy as np
2. import matplotlib.pyplot as plt
3. fig = plt.figure(1)
4. x = np.arange(10.0) / 10.0
5. y = (x + 0.1)**2
6.  
7. plt.errorbar(x, y, xerr=0.1, xlolims=True)
8. y = (x + 0.1)**3
9.  
10. plt.errorbar(x + 0.6, y, xerr=0.1, xuplims=upperlimits, xlolims=lowerlimits)
11.  
12. y = (x + 0.1)**4
13. plt.errorbar(x + 1.2, y, xerr=0.1, xuplims=True)
14.  
15. plt.xlim(-0.2, 2.4)
16. plt.ylim(-0.1, 1.3)
17. plt.savefig('Errobars_xonly.png', format="png", bbox_inches='tight')
18. plt.show()

Example 3: y errobars only

Code: [Select all] [Expand/Collapse]

1. fig = plt.figure(0)
2. x = np.arange(10.0)
3. y = np.sin(np.arange(10.0) / 20.0 * np.pi)
4.  
5. plt.errorbar(x, y, yerr=0.1)
6.  
7. y = np.sin(np.arange(10.0) / 20.0 * np.pi) + 1
8. plt.errorbar(x, y, yerr=0.1, uplims=True)
9.  
10. y = np.sin(np.arange(10.0) / 20.0 * np.pi) + 2
11.  
12. upperlimits = np.array([1, 0] * 5)
13. lowerlimits = np.array([0, 1] * 5)
14. plt.errorbar(x, y, yerr=0.1, uplims=upperlimits, lolims=lowerlimits)
15. plt.savefig('Errobars_yonly.png', format="png", bbox_inches='tight')
16. plt.xlim(-1, 10)

See refs [1], [2].

Thus, a linear regression model must now include this error $ϵ$ into it. A simple linear regression model is then given by

$$\begin{array}{}\text{(2)}& \hat{y}=b_0+b_{1}x+ϵ.\end{array}$$

This redefined linear regression model can then be used to derive an equation of a line that minimizes these errors. The approach is based on the Least Square Method. The least square method aims to minimize the errors, which are the differences between the actual values (scatter plot points) and estimated values (points along the regression line) and make them small as possible. We show how to use the ordinary least squares (OLS) method below to analytically solve for the slope and intercept estimates and use them to fit a new regression line.


Simple Linear Regression Model - Ordinary Least Squares (OLS): Closed Form Solution

There are inbuilt functions in most of the programming languages such as Python, R, MATLAB and Octave for computing the slope and the $y$-intercept estimates from the two-dimensional data points to model a linear regression line. However, we can obtain a closed form solution by defining and minimizing some objective functional. This is done by finding the estimated values ${\hat{b}}_o$ and ${\hat{b}}_1$ for the parameters $b_o$ and $b_1$ which would provide the best fit for the given two-dimensional dataset. The "best fit" in this case is interpreted as in the least-squares approach - a line that minimizes the sum of squared residuals, ${\widehat{ϵ}}_i$ (the differences between actual and predicted/estimated values of the dependent variable $y$). We can respectively, obtain the $y$-intercept and the slope estimates, ${\hat{b}}_o$ and ${\hat{b}}_1$ by taking the partial derivative of the objective/cost functional with respect to each term, setting them to zero, and calculate the slope and the $y$-intercept, where for the line of a linear regression given by the formula,
$$\begin{array}{}\text{(3)}& y_i=b_0+b_1{x}_i+{ϵ}_i,\end{array}$$

we need to find

$$\begin{array}{}\text{(4)}& {\mathrm{m}\mathrm{i}\mathrm{n}}_{{\mathrm{b}}_{\mathrm{o}},{\mathrm{b}}_1}\mathrm{F}({\mathrm{b}}_{\mathrm{o}},{\mathrm{b}}_1)\end{array}$$

where

$$\begin{array}{}\text{(5)}& F(b_o,b_1)& =\sum _{i=1}^N({\widehat{ϵ}}_i{)}^2\text{(6)}& & =\sum _{i=1}^N(y_i-b_o-b_1{x}_i{)}^2.\end{array}$$

$\sum _{i=1}^N({\widehat{ϵ}}_i{)}^2$ is called the residual sum of squares (RSS), also known as the sum of squared residuals (SSR) or the sum of squared errors (SSE) of prediction, it explains the deviations of the predicted values from actual empirical values of data. We use this as an optimality criterion for parameter and regression model selection.

We expand the cost functional $F$ to get quadratic expressions in $y$-intercept $b_o$ and slope $b_1$,

$$\begin{array}{}\text{(7)}& F(b_o,b_1)=\sum _{i=1}^N(y_i^2-2b_o{y}_i-2b_1{x}_i{y}_i+b_o^2+2b_o{b}_1{x}_i+b_1^2{x}_i^2),\end{array}$$

and then minimize it by taking the partial derivatives of the cost functional $F$ with respect to $b_o$ and $b_1$, setting the results to zero, where finally we have

$$\begin{array}{}\text{(8)}& F_{b_o}& =\sum _{i=1}^N(-2y_i+2b_o+2b_1{x}_i)=0,\text{(9)}& & ⟹\sum _{i=1}^N(y_i-b_o-b_1{x}_i)=0.\end{array}$$

Therefore,

$$\begin{array}{}\text{(10)}& {\hat{b}}_o=\overline{y}-{\hat{b}}_1\overline{x}={\displaystyle \frac{\sum _{i=1}^N{y}_i}{N}}-{\displaystyle \frac{{\hat{b}}_1\sum _{i=1}^N{x}_i}{N}}.\end{array}$$

Similarly,

$$\begin{array}{}\text{(11)}& F_{b_1}& =\sum _{i=1}^N(-2x_i{y}_i+2b_o{x}_i+2b_1{x}_i^2)=0,\text{(12)}& \sum _{i=1}^N{b}_1{x}_i^2& =\sum _{i=1}^N(y_i{x}_i-b_o{x}_i),\text{(13)}& b_1\sum _{i=1}^N{x}_i^2& =\sum _{i=1}^N{y}_i{x}_i-({\displaystyle \frac{\sum _{i=1}^N{y}_i}{N}}-b_1{\displaystyle \frac{\sum _{i=1}^N{x}_i}{N}})\sum _{i=1}^N{x}_i,\text{(14)}& b_1(\sum _{i=1}^N{x}_i^2-{\displaystyle \frac{\sum _{i=1}^N{x}_i\sum _{i=1}^N{x}_i}{N}})& =\sum _{i=1}^N{y}_i{x}_i-\sum _{i=1}^N{y}_i\sum _{i=1}^N{x}_i,\end{array}$$

which implies

$$\begin{array}{r}\text{(15)}& {\hat{b}}_1& ={\displaystyle \frac{\sum _{i=1}^N{y}_i{x}_i-{\displaystyle \frac{\sum _{i=1}^N{y}_i\sum _{i=1}^N{x}_i}{N}}}{\sum _{i=1}^N{x}_i^2-{\displaystyle \frac{\sum _{i=1}^N{x}_i\sum _{i=1}^N{x}_i}{N}}}}\text{(16)}& & ={\displaystyle \frac{\sum _{i=1}^N{y}_i{x}_i-\overline{y}N\overline{x}-\overline{x}N\overline{y}+N\overline{x}\overline{y}}{\sum _{i=1}^N{x}_i^2-2{\displaystyle \frac{\sum _{i=1}^N{x}_i\sum _{i=1}^N{x}_i}{N}}+{\displaystyle \frac{\sum _{i=1}^N{x}_i\sum _{i=1}^N{x}_i}{N}}}}\\ \text{(17)}& & ={\displaystyle \frac{\sum _{i=1}^N{y}_i{x}_i-\overline{y}\sum _{i=1}^N{x}_i-\overline{x}\sum _{i=1}^N{y}_i+N\overline{x}\overline{y}}{\sum _{i=1}^N{x}_i^2-2{\displaystyle \frac{\sum _{i=1}^N{x}_i\sum _{i=1}^N{x}_i}{N}}+{\displaystyle \frac{N\overline{x}\sum _{i=1}^N{x}_i}{N}}}}\\ \text{(18)}& & ={\displaystyle \frac{\sum _{i=1}^N{y}_i{x}_i-\overline{y}\sum _{i=1}^N{x}_i-\overline{x}\sum _{i=1}^N{y}_i+\sum _i^N\overline{x}\overline{y}}{\sum _{i=1}^N{x}_i^2-2{\displaystyle \frac{\sum _{i=1}^N{x}_i\sum _{i=1}^N{x}_i}{N}}+N{\overline{x}}^2}}\\ \text{(19)}& & ={\displaystyle \frac{\sum _{i=1}^N(x_i-\overline{x})(y_i-\overline{y})}{\sum _{i=1}^N{x}_i^2-2\overline{x}\sum _{i=1}^N{x}_i+\sum _{i=1}^N{\overline{x}}^2}}\\ \text{(20)}& & ={\displaystyle \frac{\sum _{i=1}^N(x_i-\overline{x})(y_i-\overline{y})}{\sum _{i=1}^N(x_i-\overline{x}{)}^2}}\end{array}$$

Now, we have the $y$-intercept estimate,

$$\begin{array}{}\text{(21)}& {\hat{b}}_o=\overline{y}-{\hat{b}}_1\overline{x}={\displaystyle \frac{\sum _{i=1}^N{y}_i}{N}}-{\displaystyle \frac{{\hat{b}}_1\sum _{i=1}^N{x}_i}{N}}\end{array}$$

and the gradient estimate,

$$\begin{array}{}\text{(22)}& {\hat{b}}_1={\displaystyle \frac{\sum _{i=1}^N(x_i-\overline{x})(y_i-\overline{y})}{\sum _{i=1}^N(x_i-\overline{x}{)}^2}}.\end{array}$$

Here,

• $x_i$ and $y_i$ are each measurements of two separate attributes of the same observation,
• $\overline{x}$ and $\overline{y}$ are respectively, the averages of $x_i$ and $y_i$.

$(x-\overline{x})(y-\overline{y})$ is intuitively a vector cross product, and the role of $(x_i-\overline{x})$ and $(y_i-\overline{y})$ is to center values of $x$ and $y$ at their means (averages). The multiplication and summation,

$$\begin{array}{}\text{(23)}& \sum _{i=1}^N(x_i-\overline{x})(y_i-\overline{y})\end{array}$$

is the dot product of the two vectors $(x-\overline{x})$ and $(y-\overline{y})$ that tells how parallel the two vectors are to each other. More precisely, the dot product of vectors $x$ and $y$ can be interpreted geometrically as the length of the projection of $x$ onto the unit vector $y$ when the two vectors are placed in such a way their tails coincide (see here).

Now, we're about to introduce the notion of covariance. We can view covariance as the average of the products of the mean deviations $(x_i-\overline{x})$ and $(y_i-\overline{y})$, i.e.,

$$\begin{array}{}\text{(24)}& {\displaystyle \frac{1}{N-1}}\sum _{i=1}^N(x_i-\overline{x})(y_i-\overline{y}).\end{array}$$

Covariance of two variables, $x$ and $y$ can however be defined by using the notion of expectations (weighted averages of values of $x$ and $y$,)

$$\begin{array}{}\text{(25)}& Cov(x,y)=E[(x-E(x))(y-E(y))],\end{array}$$

where $E$ is the expected value operator.

In general, covariance centers the data to the mean, and the division by $N-1$ or taking the expected values, is the process of scaling the data by the number of observations. We will see in the subsequent sections how covariance is related to correlation, which is simply a normalized version of covariance.

We can then use the Python code below to calculate the $y$-intercept estimate ${\hat{b}}_o$ and the gradient estimate ${\hat{b}}_1$.

Code: [Select all] [Expand/Collapse]

1. def simple_linear_regression(x, y):
2.     """Returns y-intercept and gradient of a simple linear regression line"""
3.     #Initial sums
4.     n = float(len(x)) #We can use x.sum()
5.     sum_x = np.sum(x)
6.     sum_y = np.sum(y)
7.     sum_xy = np.sum(x*y)
8.     sum_xx = np.sum(x**2)
9.     #Formula for gradient/slope, b1
10.     gradient = (sum_xy - (sum_x*sum_y)/n)/(sum_xx - (sum_x*sum_x)/n)
11.     #Formula for y-intercept, b0
12.     intercept = sum_y/n - gradient*(sum_x/n)
13.     return (intercept, gradient)
14.  
15. x = np.array([1,2,3,4,5])
16. y = np.array([2,4,5,4,5])
17. print("The y-intercept and gradient are:", simple_linear_regression(x, y))

The output is

Code: [Select all] [Expand/Collapse]

1. The y-intercept and gradient are: (2.2, 0.6)

Initially, we plotted the regression line by automatically obtaining the slope and the intercept from a call to python function, numpy.polyfit. Since we now know the slope and intercept of our regression line fitting function, we can use them to fit the regression line, as shown in the code below.

Code: [Select all] [Expand/Collapse]

1. #Import required libraries
2. import numpy as np
3. import matplotlib.pyplot as plt
4. #Create data points
5. x = np.array([1,2,3,4,5])
6. y = np.array([2,4,5,4,5])
7.      
8. fig  = plt.subplot(111, aspect='auto')
9. plt.title("Regression line: b0 and b1 obtained analytically")
10. plt.xlabel("Independent variable x")
11. plt.ylabel("Dependent variable  y")
12. plt.ylim(top=7)
13. fig.plot(x, y, 'o', markersize=7, color='black', alpha=1.0, label="Scatter plot")
14. plt.axvline(x=np.mean(x), color='red', linestyle='--', alpha=1.0, label="x Mean")
15. plt.axhline(y=np.mean(y), color='green', linestyle='--', alpha=1.0, label="y Mean")
16. #Define the slope and intercept of the fit function
17. b1, b0 = 0.6, 2.2
18. f = lambda x: b1*x + b0
19. plt.plot(x,f(x), c="purple", linestyle = "-", label="Regression line")
20. #Plot half errorbars
21. plt.errorbar(x,f(x),xerr=0, yerr=[np.zeros(5), y-f(x)], fmt="o", c="blue", label=r"${\rm Error,} \ \epsilon$")
22. plt.legend(loc=2)
23. plt.savefig('new_fit_regression_line.png', format="png", bbox_inches='tight')

The result is the same as before:



Relationship Between Covariance and Correlation

From the gradient estimate we obtained above,

$$\begin{array}{}\text{(26)}& {\hat{b}}_1& ={\displaystyle \frac{\sum _{i=1}^N(x_i-\overline{x})(y_i-\overline{y})}{\sum _{i=1}^N(x_i-\overline{x}{)}^2}}\text{(27)}& & ={\displaystyle \frac{\frac{1}{N-1}\sum _{i=1}^N(x_i-\overline{x})(y_i-\overline{y})}{\frac{1}{N-1}\sum _{i=1}^N(x_i-\overline{x}{)}^2}}\text{(28)}& & ={\displaystyle \frac{Cov(x,y)}{{\sigma }^2(x)}}\text{(29)}& & ={\displaystyle \frac{Cov(x,y)}{Var(x)}}\text{(30)}& & =\frac{\mathrm{S}\mathrm{a}\mathrm{m}\mathrm{p}\mathrm{l}\mathrm{e}\mathrm{c}\mathrm{o}\mathrm{v}\mathrm{a}\mathrm{r}\mathrm{i}\mathrm{a}\mathrm{n}\mathrm{c}\mathrm{e}\mathrm{o}\mathrm{f}x\text{and}y}{\mathrm{S}\mathrm{a}\mathrm{m}\mathrm{p}\mathrm{l}\mathrm{e}\mathrm{v}\mathrm{a}\mathrm{r}\mathrm{i}\mathrm{a}\mathrm{n}\mathrm{c}\mathrm{e}\mathrm{o}\mathrm{f}x}.\end{array}$$

We know that
$$\begin{array}{}\text{(31)}& \mathrm{c}\mathrm{o}\mathrm{r}\mathrm{r}\mathrm{e}\mathrm{l}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n},r_{xy}& ={\displaystyle \frac{Cov(x,y)}{\sqrt{Var(x)Var(y)}}}={\displaystyle \frac{Cov(x,y)}{\sqrt{Var(x)}\sqrt{Var(y)}}}\text{(32)}& & =\sqrt{{\displaystyle \frac{Var(x)}{Var(y)}}}{\displaystyle \frac{Cov(x,y)}{Var(x)}}.\end{array}$$

$$\begin{array}{}\text{(33)}& \mathrm{T}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{e},{\hat{b}}_1={\displaystyle \frac{Cov(x,y)}{Var(x)}}=r_{xy}\sqrt{{\displaystyle \frac{Var(y)}{Var(x)}}},\text{(34)}& {\hat{b}}_1=r_{xy}{\displaystyle \frac{s_y}{s_x}}.\end{array}$$

We have introduced new parameters as follows:

• $Cov(x,y)$, sample covariance of $x$ and $y$.
• $Var(x)$ and $Var(y)$ are sample variances of $x$ and $y$. If a set of data points has many values far away from its mean, it will have a high variance and vice versa.
• $r_{xy}$, the sample correlation coefficient between $x$ and $y$.
• $s_x$ and $s_y$, the sample standard deviations of $x$ and $y$. The standard deviation is a measure that is used to quantify the amount/degree of variation or dispersion (how spread) of the set of data values. A small standard deviation means that the data points tend to be close to (clustered around) the mean of the set, while a larger standard deviation means the data points are spread far from the mean.

Furthermore,

• The variance is a covariance of a variable, say $x$, with itself, i.e. $Cov(x,x)=Var(x)$.
• The square root of variance is called standard deviation.

Covariances are very useful in many ways, but, the problem with covariances is that they are hard to use for comparison reasons, for example, if we want to use them to establish how two quantities, such as a set of weights (in kilograms) and a set of heights (in metres) covary. In this case, you are subject to get different covariances when you repeat calculations using the same observational data or measurements but this time in different units, e.g., grams and centimetres, instead of previous units of kilograms and metres. It will however be hard to tell, for example, if height and weight covariate better than, e.g. the length of your arms and legs, simply because of different scales used to calculate the covariances. Due to these reasons, covariance is unbounded and lacks a proper context when comparing to other covariances.

The solution to these challenges is to normalize the covariance, that's dividing the covariance, $Cov(x,y)$ by the square root of the product of the variances (standard deviations) of the variables/vectors $x$ and $y$. This adjustment represents the diversity and scale in both the covariates, resulting in a value that is certainly between -1 and 1, namely; the correlation coefficient, $r_{xy}$. Whatever unit the original variables were in, we are assured to always get the same result, and this will also ensure that we can, to a certain extent, compare whether two variables correlate more than two others, simply by comparing their correlation. Very often, correlation, obtained by normalizing (adjusting or standardizing) covariance, is used to compare datasets more easily.

We have seen that covariance simply centers the data. Correlation (dimensionless), however, not only centers the data but also normalizes the data using the standard deviations of the respective variables. Here are few remarks we can make in regard to correlation:

• Positive correlation/relationship means the value of dependent variable increases/decreases as the value of independent variable increases/decreases ($+1$ means perfectly positively correlated).
• Negative correlation/relationship means the value of dependent variable decreases/increases as the value of independent variable increases/decreases ($-1$ means perfectly inversely correlated).
• If we get the value of correlation close to zero, it means that the two variables are uncorrelated (zero means perfectly uncorrelated).

The relationship between the Correlation Coefficient, $r_{xy}$ and the Coefficient of Determination, $R^2$

From
$$\begin{array}{}\text{(35)}& \hat{y}={\hat{b}}_o+{\hat{b}}_{1}x\end{array}$$,

Let's substitute the expressions we obtained for the estimates of $b_o$ and $b_1$ into the above equation to have,

$$\begin{array}{}\text{(36)}& \hat{y}& =\overline{y}-{\hat{b}}_1\overline{x}+r_{xy}{\displaystyle \frac{s_y}{s_x}}x=\overline{y}-r_{xy}{\displaystyle \frac{s_y}{s_x}}\overline{x}+r_{xy}{\displaystyle \frac{s_y}{s_x}}x,\text{(37)}& \hat{y}-\overline{y}& =r_{xy}{s}_y{\displaystyle \frac{x-\overline{x}}{s_x}},\text{(38)}& ⟹{\displaystyle \frac{\hat{y}-\overline{y}}{s_y}}& =r_{xy}{\displaystyle \frac{x-\overline{x}}{s_x}}.\end{array}$$

This last equation is in the form $y=b_o+b_{1}x$ (linear equation), where $r_{xy}$ is the slope of the regression line of the standardized data points. Here, the $y$-intercept is zero, and hence this regression line passes through the origin.

We can now generalize the notion, $\overline{x}$, to indicate the average value over the set of random samples, $x$ and $y$ as

$$\begin{array}{}\text{(39)}& \overline{xy}={\displaystyle \frac{1}{N}}\sum _{i=1}^N{x}_i{y}_i.\end{array}$$

This notion allows us to construct a formula for $r_{xy}$, that's

$$\begin{array}{}\text{(40)}& r_{xy}={\displaystyle \frac{\overline{xy}-\overline{x}\overline{y}}{\sqrt{(\overline{x^2}-{\overline{x}}^2)(\overline{y^2}-{\overline{y}}^2)}}}\end{array}$$.

Therefore, we can deduce that, the coefficient of determination, $R^2$ is equal to the correlation coefficient squared, $r_{xy}^2$ for the simple linear regression model with a single independent variable. How? Let's verify this claim mathematically.

From

$$\begin{array}{}\text{(41)}& {\displaystyle \frac{\hat{y}-\overline{y}}{s_y}}& =r_{xy}{\displaystyle \frac{x-\overline{x}}{s_x}},\text{(42)}& r_{xy}& ={\displaystyle \frac{s_x}{s_y}}\times {\displaystyle \frac{\hat{y}-\overline{y}}{x-\overline{x}}},\text{(43)}& r_{xy}^2& ={\displaystyle \frac{s_x^2}{s_y^2}}{\displaystyle \frac{(\hat{y}-\overline{y}{)}^2}{(x-\overline{x}{)}^2}},\text{(44)}& r_{xy}^2& ={\displaystyle \frac{{\displaystyle \frac{\sum _{i=1}^N(x_i-\overline{x}{)}^2}{N-1}}}{{\displaystyle \frac{\sum _{i=1}^N(y_i-\overline{y}{)}^2}{N-1}}}}\times {\displaystyle \frac{(\hat{y}-\overline{y}{)}^2}{(x-\overline{x}{)}^2}}.\end{array}$$

Introducing indices and summation for $(\hat{y}-\overline{y}{)}^2$ and $(x-\overline{x}{)}^2$, we have

$$\begin{array}{}\text{(45)}& ⟹r_{xy}^2& ={\displaystyle \frac{{\displaystyle \frac{\sum _{i=1}^N(x_i-\overline{x}{)}^2}{N-1}}}{{\displaystyle \frac{\sum _{i=1}^N(y_i-\overline{y}{)}^2}{N-1}}}}\times {\displaystyle \frac{\sum _{i=1}^N({\hat{y}}_i-\overline{y}{)}^2}{\sum _{i=1}^N(x_i-\overline{x}{)}^2}},\text{(46)}& r_{xy}^2& ={\displaystyle \frac{\sum _{i=1}^N({\hat{y}}_i-\overline{\hat{y}}{)}^2}{\sum _{i=1}^N(y_i-\overline{y}{)}^2}}=r_{y\hat{y}}^2=R^2.\end{array}$$

This analysis can be extended to higher dimensions than 2.

Computing $R^2$

In regression analysis, $R^2$ tells us how well a regression line predicts the estimates of the actual values. It compares two distances: the distances of actual values from the mean to the distances of estimated values from the mean of the dependent variable $y$. To calculate $R^2$ from the given dataset, we need to determine:

• the mean of actual values of the dependent variable; in this case, we will use values for the variable $y$ from the table above with the means, $\overline{y}$ of $4$;
• distances of actual values for the depependent variable, $y$ to the mean, $y-\overline{y}$;
• distances of estimated values (regression line) to the mean, $\hat{y}-\overline{\hat{y}}$.

Finally, we compare distances of the actual values to the mean to the the distances of the estimates to the mean, i.e., $(y-\overline{y}{)}^2$ and $(\hat{y}-\overline{\hat{y}}{)}^2$ using the expression,

$$\begin{array}{}\text{(47)}& R^2={\displaystyle \frac{(\hat{y}-\overline{\hat{y}}{)}^2}{(y-\overline{y}{)}^2}}.\end{array}$$

We summarize this process using the table below.

$$\begin{array}{|ccccccc|}\hline x_i& y_i& y_i-{\overline{y}}_i& (y_i-{\overline{y}}_i{)}^2& {\hat{y}}_i& {\hat{y}}_i-{\overline{\hat{y}}}_i& ({\hat{y}}_i-{\overline{\hat{y}}}_i{)}^2\\ 1& 2& -2& 4& 2.8& -1.2& 1.44\\ 2& 4& 0& 0& 3.4& -0.6& 0.36\\ 3& 5& 1& 1& 4.0& 0.0& 0.0\\ 4& 4& 0& 0& 4.6& 0.6& 0.36\\ 5& 5& 1& 1& 5.2& 1.2& 1.44\\ \hline\end{array}$$

Note that, the estimated values, $\hat{y}$, are computed using the equation $\hat{y}=b_o+b_{1}x=2.2+0.6x$, where $b_0=2.2$ and $b_1=0.6$ as we have previously seen. $\hat{y}$, a predicted value of $y$ for a given value of $x$, is the line that minimizes the sum of squared errors, ${ϵ}^2$.

We can therefore compute R-squared as follows:

$$\begin{array}{}\text{(48)}& R^2& ={\displaystyle \frac{(\hat{y}-\overline{\hat{y}}{)}^2}{(y-\overline{y}{)}^2}}={\displaystyle \frac{\sum _{i=1}^N({\hat{y}}_i-\overline{\hat{y}}{)}^2}{\sum _{i=1}^N(y_i-\overline{y}{)}^2}}\text{(49)}& & ={\displaystyle \frac{1.44+0+0.36+1.44}{4+0+1+0+1}}={\displaystyle \frac{3.6}{6}}=0.6.\end{array}$$

$R^2$ is the proportion of the variation in $y$ being explained by the variation in $x$.

In most cases, the formula for R-squared is commonly acronymised as

$$\begin{array}{}\text{(50)}& R^2& ={\displaystyle \frac{\mathrm{E}\mathrm{S}\mathrm{S}/\mathrm{S}\mathrm{S}\mathrm{R}}{\mathrm{T}\mathrm{S}\mathrm{S}/\mathrm{S}\mathrm{S}\mathrm{T}}}\text{(51)}& & ={\displaystyle \frac{\sum _{i=1}^N({\hat{y}}_i-\overline{\hat{y}}{)}^2}{\sum _{i=1}^N(y_i-\overline{y}{)}^2}},\end{array}$$

where,

• ESS/SSR means ESS or SSR, similarly, TSS/SST means TSS or SST.
• ESS is the explained sum of squares, also known as the model sum of squares or sum of squares due to regression ("SSR" – not to be confused with the sum of squared residuals (SSR) that we will explain below).
• The residual sum of squares (RSS), $\sum _{i=1}^N(y_i-{\hat{y}}_i{)}^2$,
  
  also known as the sum of squared residuals (SSR) ( not to be confused with the sum of squares due to regression "SSR" above) or the sum of squared errors of prediction (SSE). This is the sum of the squares of residuals (sum of deviations of predicted values/estimates from actual empirical/observational data).
• TSS is the total sum of squares, alternatively known as sum of the squares (SST), defined as the sum, over all observations, of the squared differences of each observation from the overall mean.
• Generally (see here), the total sum of squares (TSS) = explained sum of squares (ESS) + residual sum of squares (RSS).

Most importantly, don't care much about these terms and their abbreviations, focus on mathematics and how we arrived to $R^2$.

Enjoy maths

Welcome for suggestions and improvements: @Simulink @cache @jiwe @FantasmaDiamato @kipanga @RealityKing @Forbidden_Technology

[RealityKing]

Thank you very much @Eli , great post!

[Joseph Bundala]

Well, this is a very good topic @Eli and well presented. In Control Systems, it is called System Identification. However, in industrial processes these algorithms are used to find coefficients of very higher order/exponential/discrete-time variant non-linear data from using only input (u) and output (y) data from a system. I am introducing another algorithm called Recursive Least Square (RLS) which estimates the coefficients of a polynomial recursively. Unlike LSM as described by @Eli which just computes once the values of coefficients. RLS enables to compute error divergences towards the best solution while updating coefficients.

Assume a Telecommunication channel which uses TDM E1, voice through it will have noise, disturbances and floating voltages accompanied to the receiver. RLS algorithm is used to perform noise cancellation in GSM networks by smoothing the signals through Signal Processing.
I will demonstrate this example of data which has few iterations.

u	y
uu1	yy1.9
uu2	yy4
uu3	yy6.2

Algorithm

Consider an input signal ,output signal whole relationship is linear such that





Then, the coefficient of a linear system can be obtained by:-



Initialization of values;

, this is the best known value of
, diagonal matrix,

Plots for different values of , and MSE . Is MSE expected to decrease or increase towards the best solution?

RLS and LSM

MSE

Matlab Code

Code: [Select all] [Expand/Collapse]

1. function RLS
2.  
3.   clc
4.   clear all
5.   close all
6.  
7.   u = [ 1, 2, 3]';
8.   y = [1.9, 4, 6.2]';
9.   %u = [1, 2, 3, 4,5]';
10.   %y = [2, 4, 5, 4, 5]';
11.  
12.   %x = ones(length(u),1);
13.   b(1) = 0;     % initial value of b
14.   p(1) = 100;  
15.   lamda = 0.98;
16.   %RLS- Recursive Least Square
17.  
18.   scatter(u,y)
19.   hold on
20.   %uu = [x u];
21.   %bb = inv(uu'*uu)*uu'*y;
22.   bb = inv(u'*u)*u'*y
23.   plot(u,u*bb,'-p')
24.  
25.   error = y;
26.   %u = [x,u];
27.   for j = 1: length(y)
28.     k  = p(j)*u(j)/(1+u'(j)*p(j)*u(j));
29.     p(j+1)  = (p(j) - ((p(j)*u(j)*u'(j)*p(j))/(lamda+(u'(j)*p(j)*u(j)))))/lamda;
30.     b(j+1) = b(j) + (k* (y(j) - (u'(j)*b(j))));
31.     error(j) = y(j)-b(j)*u(j);
32.     Z = [u y b(j+1)*u error];
33.     z = dataframe (Z, "colnames", {"u", "y", "y_estimated","error"})
34.     MSE(j) = 1/length(u)* sum( error.^2);  
35.   endfor
36.   grid on  
37.   plot( u, u*b(2),'--',u, u*b(3),u,u*b(4))
38.   text(2.1253,4.0036,'\leftarrow b = 1.8812')
39.   text(2.3097,4.5509,'\leftarrow b = 1.9761')
40.   text(2.6138,5.3277,'\leftarrow b = 2.0343')
41.   legend("Scatter","LSM","RLS First Estimation","RLS Second Estimation","RLS Third Estimation","location",'north')
42.   h=ylabel("y")
43.   xlabel("u")
44.   title("y = bu")
45.   grid on
46.   set(get(gca,'ylabel'),'rotation',0)
47.   saveas(1,"RLS.png")
48.  
49.   figure
50.   subplot(2,1,1)
51.   semilogy(abs(error))
52.   ylabel("error")
53.   title("log error")
54.   set(get(gca,'ylabel'),'rotation',0)
55.   grid on
56.  
57.   subplot(2,1,2)  
58.   plot(MSE)
59.   title("Mean Square Error(MSE)")
60.   ylabel("MSE")
61.   set(get(gca,'ylabel'),'rotation',0)
62.   grid on
63.   saveas(2,"MSE.png")
64.    
65. endfunction
66.  

To estimate a function with data points given by Eli is a bit tricky, these initializations will help you estimate those coefficients by RLS:

Code: [Select all] [Expand/Collapse]

1. %u = input signal
2. %M =  number of coefficients
3. %N = length of data
4. b = zeros(M,1);
5. p = 100*cov(u).^2*eye(2);
6. for j = M:N
7.     newMat = u(j:-1:j-M+1);
8.     .
9.     .
10.     k = .....
11.     p = .....
12.     b = ....
13. endfor
14. plot blah blah
15. }

References

S. Haykin- Adaptive Filtering Theory - Prentice Hall, 2002.

[Eli]

Thank you @Simulink for your contribution. I learned something interesting, a simple idea from which linear regression arises can be amplified to do pretty cool things!

[FantasmaDiamato]

I translate it artistically......and woo it's mind blowing

[Eli]

I translate it artistically......and woo it's mind blowing

A person who is doing linear regression analysis can seamlessly apply these concepts and get what he/she wants out of the box.

Nice seeing you here @FantasmaDiamato

[Eli]

Here are nice references for Regression Analysis:

1. http://cds.cern.ch/record/302220/files/9605002.pdf

2. https://pdfs.semanticscholar.org/7a07/5 ... 0f687e.pdf

[Eli]

Here is a related topic

Data Representation and Visualization: From Scatter Plot to 2-Dimensional Histogram to Contours

Appropriate representation of data is a key point in data analysis. Hundreds to thousands of numerical data points, for example, can be difficult to interpret unless correctly visualized. Improper choice of visualization can however lead to overplotting a

[Eli]

If we consider a regression model with one independent variable, $y={\beta }_{1}x+{\beta }_2$ (see post #1), the coefficient of determination $R^2$ is equivalent to the Squared Pearson Correlation Coefficient between the observed values $y_i$ and the fitted/estimated values ${\hat{y}}_i$.

[Eli]