I am stuck trying to understand the basic calculation of ordinary least squares. From Wikipedia:
$$y = \beta X^T + \varepsilon$$
where $X$ is the independent variable, $Y$ is the dependent variable and $X^T$ denotes the transpose of $X$.
Why are we taking the transpose? In a simple linear function $f(x,y):y= ax + b$, we do not use the transpose of $x$, we just use the variable $X$ as is. So I am not sure why in linear regression we need to take the transpose of $X$?
Let me begin by saying it seems as though you're not being very careful when you write mathematical expressions. Nevertheless, I think I can understand what you're asking, so I will try to answer.
To begin, we ought to straighten out the notation being used. Doing this will enable us to communicate more clearly and identify any misunderstandings that may exist.
When you write $f(x,y):y= ax + b$, I take it that you really mean the following $$ y = ax + b $$ where $a$ and $b$ are constants and real numbers (although, strictly speaking they need not be real numbers). Furthermore, $y$ is a function of only one variable, $x$, so we write $f(x)$ instead of $f(x,y)$. Pay particular attention to the fact that the commutative property of multiplication is satisfied when multiplying both $a$ and $x$ together. In other words, the order in which $a$ and $x$ are multiplied does not matter.
Now enter the other equation; namely, what you've written as $y = \beta X^{T} + \epsilon$. This is actually matrix notation and it has a different meaning to the (more everyday?) algebraic notation shown in the previous paragraph. The linear regression model can, in fact, be more commonly written in the following form $$ y = X \beta + \epsilon $$ where $y$ is a $(T \times 1)$ vector, $X$ is a $(T \times K)$ matrix, $\beta$ is a $(K \times 1)$ vector, and $\epsilon$ is a $(T \times 1)$ vector. The notation being used in brackets refers to the dimensions of each vector or matrix, where the first number refers to the number of rows and the second to the number of columns. So, for example, $X$ is a matrix with $T$ rows and $K$ columns.
Note the distinction! $X$ and $\beta$ are matrices not like $a$ and $x$, which are just scalar values. Crucially, the commutative property of (matrix) multiplication is not always satisfied in the world of linear algebra; the general case is that it is not satisfied. In other words, in general, $X\beta \neq \beta X$. Furthermore, such operation may not even be permitted, since in matrix algebra, matrices must be conformable for multiplication if they are to be multiplied.
In order for two matrices to be conformable for matrix multiplication, we say that the number of columns of the left matrix must be the same as the number of rows of the right matrix. In the linear regression model, $X \beta$ is possible because $X$, the left matrix, has $K$ columns and $\beta$, the right matrix, has $K$ rows. On the other hand, $\beta X$ would not be possible because $\beta$, the first matrix, has $1$ column while $X$, the second matrix, has $T$ rows - unless, of course, $T = 1$.
At this stage, please revise what you meant when you wrote: $y = \beta X^{T} + \epsilon$.
Lastly, the Wikipedia page that you link us to shows that the linear regression model can also be written in a form involving a transpose: $$ y_{i} = x_{i}^{T} \beta + \epsilon_{i} $$ where both $x_{i}$ and $\beta$ are column vectors of dimension $(p \times 1)$.
If anything I've written has made sense, you should be able to decipher why $x_{i}^{T} \beta$ and not $x_{i} \beta$.
Graeme Walsh has an excellent answer here. However, I decided I would take a crack at it from a combination of linear algebra and programming in R to show why this matters (and so perhaps I can use this to reference people to down the road). Before reading this, if you don't know any linear algebra, this link may help you follow along. I'm guessing you already know a lot of this, but for clarity I explain everything to its basics (at least as much as writing this answer won't be too long).
You may be familiar with the lm function in R, which estimates a linear model by default in the OLS way. I will fit this linear regression in R first from the well-known iris data, which has measures of flower dimensions such as the widths and lengths of petals and sepals. We can specify a simple model here with one intercept and two slopes as follows:
$$ \text{Petal Length} = \beta_0 + {\beta_1} \text{Petal Width} + \beta_2 \text{Sepal Width} + \epsilon $$
First the model fit, which we save as an object called fit:
fit <- lm(Petal.Length ~ Petal.Width + Sepal.Width,
iris)
We won't peek at this just yet. We will come back to this at the end, but we have saved this regression output here for now to compare against our hand-calculated version later. Now on to the explaining...as Graham already alluded to, we do not estimate $y$ by transposing the vector $x$, as all this does is flip the vector:
$$ \text{From this} = \begin{bmatrix} x_1 \\ x_2 \\ . \\ . \\ . \\ x_i \end{bmatrix} % \hspace{5mm} \text{To this} = \begin{bmatrix} x_1 ,x_2, . . . x_i \end{bmatrix} $$
Multiple regression can potentially have a large number of predictors, so it is more efficient to use matrices to define the regression instead. This is also because we have to solve a set of linear equations such as below, where each $x_i$ is a value in the vector $x$ (and this can of course be generalized to more predictors):
$$ y = \beta_0 + \beta_1 x_1 + \epsilon \\ y = \beta_0 + \beta_1 x_2 + \epsilon \\ .\\ .\\ .\\ y = \beta_0 + \beta_1 x_i + \epsilon $$
To generalize to multiple linear regression, we need to first consider a matrix which combines our intercept with our raw data from each predictor called $X$ (note the capitalization here signifies in this context the matrix). Our $X$ in our R case contains three elements for our regression:
This can be easily created in R like so, where cbind simply glues these three columns together (the $1$ will simply repeat here $150$ times to match the other two columns):
X <- cbind(
1,
iris$Petal.Width,
iris$Sepal.Width
)
Printing X gets us this matrix (I show the first ten rows here to save space):
[,1] [,2] [,3]
[1,] 1 0.2 3.5
[2,] 1 0.2 3.0
[3,] 1 0.2 3.2
[4,] 1 0.2 3.1
[5,] 1 0.2 3.6
[6,] 1 0.4 3.9
[7,] 1 0.3 3.4
[8,] 1 0.2 3.4
[9,] 1 0.2 2.9
[10,] 1 0.1 3.1
Now we just need our $Y$, which will just be a vector of petal lengths.
Y <- iris$Petal.Length
So now we just need to get $\beta$, the vector of coefficients, by solving with the formula $y = X \beta + \epsilon$ as noted by Graeme's answer. To achieve this, we could maybe just try to multiply $X$ by $Y$ using the %*% operator for matrix multiplication:
X %*% Y
This doesn't really make much sense to do, but illustrates something about matrix multiplication, as it predictably fails:
Error in X %*% Y : non-conformable arguments
Why does this operation fail? In order to do matrix multiplication, we need the following condition: the columns of the first matrix, which we can call matrix $A$, must be equal to the rows in the second matrix, which we may call matrix $B$. So if matrix $A$ is a $2 \times 3$ matrix and matrix $B$ is a $3 \times 2$ matrix, this will work! But if $A$ was instead a $2 \times 4$ matrix, it would fail to give us a solution. How this is achieved involves creating inner products to form a solved matrix $C$, which is kindly described in the link above (along with matrix addition). As an example of that, here is the product of two matrices with the matching dimensions necessary for a solution:
$$ \begin{bmatrix} 5 & 2 \\ 3 & 8 \\ 0 & 7 \end{bmatrix} \times \begin{bmatrix} 9 & 1 & 0 \\ 0 & 5 & 4 \end{bmatrix} = \begin{bmatrix} 45 & 15 & 8 \\ 27 & 43 & 32 \\ 0 & 35 & 28 \end{bmatrix} $$
Now on to transposes and inverses. As noted earlier, a transpose of a vector simply flips a row vector into a column vector and vice versa. For the transpose of a matrix, $A^T$, we interchange the columns and rows, such as shown below:
$$ A = \begin{bmatrix} 5 & 2 \\ 3 & 8 \\ 0 & 7 \end{bmatrix} % A^T= \begin{bmatrix} 5 & 3 & 0 \\ 2 & 8 & 7 \end{bmatrix} $$
We can test this on our $X$ matrix with the following code:
X # original matrix
t(X) # transposed matrix
You'll now notice $X$ is a $3 \times 150$ matrix when it was originally a $150 \times 3$ matrix. Now the nice part about a matrix and it's transpose is this: multiplying them together gives you a symmetric matrix, but the order matters here for linear regression. If we wrap the product of $X$ and $X^T$ with the dim function, we can see the row by column dimensions of the matrix, where one becomes a $3 \times 3$ matrix while another becomes a $150 \times 150$ matrix:
dim(X %*% t(X)) # what we dont want
dim(t(X) %*% X) # what we want
We want the $3 \times 3$ solution, not the first solution. This part will be crucial later, but now we have to figure out the inverse. Remember that an inverse of a matrix gives us an identity matrix of the original, which if multiplied by the original matrix, will equal itself. This is the equivalent of "multiplying by 1" in math. Now there are different ways of achieving this, such as Gaussian or Gauss-Jordan Elimination, QR decomposition, and perhaps some other methods, but because this is a bit time consuming to write out, explain, and code, I will simply rest on the solve function to do that for me, which will calculate a solution for $X^T X$. EDIT: I show how to do this with a custom function at the end of this answer to further elucidate what is going on under the hood.
solve(t(X) %*% X)
This gives back a $3 \times 3$ matrix solution which we can now use to find out our coefficients in the linear regression:
[,1] [,2] [,3]
[1,] 0.46981785 -0.042111063 -0.134969247
[2,] -0.04211106 0.013339489 0.008540963
[3,] -0.13496925 0.008540963 0.040795612
The formula for achieving this is for our case of one intercept and two slopes would be:
$$ \beta = \begin{bmatrix} \beta_0 \\ \beta_1 \\ \beta_2 \end{bmatrix} = (X^T X)^{-1} X^T Y $$
which can be verbalized as "the inverse of the product $X$ transpose $X$ times transpose $X$ times $Y$." Now think back for a second. Our inverse of $X^T X$ was a matrix with $3$ columns, and we already know that $X^T$ has three rows, which means that we can perform matrix multiplication, something we wouldn't have been able to do if the order was flipped when finding the inverse. See below for proof:
solve(X %*% t(X)) %*% t(X) # doesn't work
solve(t(X) %*% X) %*% t(X) # works like magic
Now our final step is to simply multiply all of that by $Y$. Remember that $Y$ by itself is just a column vector of $150$ rows. But since the left hand side of the solution of $(X^T X)^{-1}$ has $150$ columns, we can again get our product! Our complete equation for solving for $\beta$ is simply this code:
betas <- solve(t(X) %*% X) %*% t(X) %*% Y
That was a lot of work, but now we can finally check to see if the regression from our saved lm object matches what we get here. As a reminder, we fit our original regression with this code:
fit <- lm(Petal.Length ~ Petal.Width + Sepal.Width,
iris)
We can now print both our estimated betas and those from the saved model using betas for the former and coef(fit) for the latter.
betas
coef(fit)
And they match!
> betas
[,1]
[1,] 2.2581635
[2,] 2.1556105
[3,] -0.3550346
> coef(fit)
(Intercept) Petal.Width Sepal.Width
2.2581635 2.1556105 -0.3550346
It may be illustrative to see how one can compute the inverse of $X^T X$ without using a simple wrapper function like solve. We can create a custom function that uses Gauss-Jordan elimination, which is annotated and explained below:
#### Gauss Jordan Inverse Function ####
gauss.jor <- function(a) {
#### Get Dimensions ####
n <- dim(a)[1] # Get number of rows/columns (it's a square matrix)
#### Augment With the Identity Matrix ####
a <- cbind(a, diag(n))
#### Loop Over Each Row of Matrix ####
for (k in 1:n) {
#### For Each Row, Loop Over All Other Rows (Our Operations) ####
for (i in 1:n) {
### Don't Change the Row We're Currently Working On ###
if (i != k) {
#### Factor Method ####
# Calculate the factor by which to multiply the current row,
# which is the element at the intersection of current row & column k,
# divided by the diagonal element of the current column.
factor <- a[i, k] / a[k, k]
#### Operations ####
# Subtract current row (multiplied by the factor) from all other rows
# This creates zeros below/above the diagonal element of current column.
a[i, ] <- a[i, ] - factor * a[k, ]
}
}
#### Morph to 1 ####
# Divide current row by its diagonal element to make the diagonal element 1,
# which is necessary for completing Gauss-Jordan elimination.
a[k, ] <- a[k, ] / a[k, k]
}
#### Return Inverse of Augmented Matrix ###
# This is the right half after the original matrix has been transformed into
# the identity matrix. If we print a we get a bunch of unnecessary 1's.
return(a[, -(1:n)])
}
#### Use the Gauss-Jordan Method ####
gauss.jor(t(X) %*% X) # original computation
gauss.jor(t(X) %*% X) %*% t(X) %*% Y # the betas
And you can see the beta coefficients are the same. Note that this method is not ideal for complex computation and is just for illustration since it is simpler than the actual solve function (and uses a different method than Gauss-Jordan). In any case, it is close enough to show how one can hand-calculate all of this.
