Least squares estimation

How do linear models estimate parameters?

In our Linear models lecture, we used the lm() function in R, but we haven't talked about how it estimates the parameters. It uses this least squares estimation, which we will now discuss the theory for.

This is the equation we use in the linear model to determine the expected value of $y$ :

E [y_{i} | x_{i}] = x_{i} β = β_{0} + x_{i 1} β_{1} + \dots

For each value of $β$ (intercept or slope for each covariate), we look at the squared residuals to evaluate how accurate it is (how much our observed value differs from the expectation):

r_{i}^{2} = (y_{i} - E [y_{i} | x_{i}])^{2} = (y_{i} - x_{i} β)^{2}

We are looking for the value of $β$ that gives us the lowest squared residuals in all observations. The easiest way to do this is to get the sum of all the squared residuals, and figure out which value of beta gives the lowest sum. We call this the optimal parameter estimate $\hat{β}$ , and its formula is:

\hat{β} = min_{β} \sum_{i = 1}^{n} (y_{i} - x_{i} β)^{2}

In the lecture notes, they give some instructions on how to manually calculate this in R. This is what's going on "under the hood" in the lm() function.

The gist of it is that one sum of squares value is given by sum(y - beta*x)^2. Using a for loop over many beta values between 0 and 2 can be used to get the minimum.

Least squares estimate vs maximum likelihood estimate

Under the Gaussian Observation Model, what they call a "maximum likelihood estimate" is the same as these least squares estimate. The proof for that is here, but I think we don't exactly need to know it.

Proof

Since it's a Gaussian observation model, we assume that the values of the response variable $y_{i}$ are normally distributed around a mean made of $x_{i} β$ , or the covariate $x$ times a parameter $β$ :

y_{i} \sim N (x_{i} β, σ^{2})

If we assume that all observations are independent, the joint probability density for all observations is given by this function with a big Pi ( $\prod$ ) product (like $\sum$ but with multiplying instead of adding.)

p (y | x) = \prod_{i = 1}^{n} N (y_{i} | x_{i} β, σ^{2})

This means that the total probability distribution is the product of the distributions of all the $y_{i}$ observations.

This gives us an idea about the uncertainty of each value of $y$ based on the covariates $x$ (since it's like a likelihood of getting each value of $y$ ). We are considering $x_{i}$ and $y_{i}$ as fixed, then choosing values of $β$ and $σ$ to maximize the joint probability density.

The maximum likelihood estimate for $β$ e.g. $\hat{β}$ is given by:

\hat{β} = max_{β} \prod_{i = 1}^{n} N (y_{i} | x_{i} β, σ^{2})

Since the normal distribution can be modelled by the probability function:

p (y_{i}) = \frac{1}{\sqrt{2 π σ^{2}}} e^{- (y_{i} - μ)^{2} / 2 σ^{2}},

Then this gets transformed into the least squares estimate from earlier. I don't really understand the math but I think it's not essential.

Note that this is only true for the Gaussian observation model - it's not true for the distribution of a t-test for example, or for other kinds of observation models used in Generalized linear models.