Probability of sample point given a Linear Regression

Question

This question may be ill-posed, but hopefully you all can help talk me through it.

Given a probability density function $f(\cdot)$ parameterized by one or more parameters $\theta$, we can compute the likelihood of a sample $x_0$ as $L(\theta ; x_0) = f(x_0 \mid \theta)$

Now, suppose that we are given a set of input data points $\{(x_1,y_1), (x_2,y_2), \ldots, (x_n,y_n)\}$ and we fit a linear regression $\hat h(x) = \hat\alpha + \hat\beta x$. Now I know that given some input $x_0$, if the assumptions of the model are satisfied, then $\hat h(x_0) = \operatorname{E}[y\mid x_0]$. This presence of an expectation makes me think there there must be some underlying probability density function $f(y\mid x)$ that I can use to compute likelihoods of arbitrary input points $(x_i,y_i)$ given the fitted regression model.

So my question is, if this probability distribution does exist, what form does it take (and how might I obtain it in R/Python), or if it doesn't exist, then where does my misunderstanding lie?

The background to all of this is that I'm trying to perform model selection between two nested linear models using a likelihood ratio test, but can't figure out how to compute the likelihoods needed to perform the test in the first place.

Answer 1

Remember that the likelihood function is the a function over the unknown parameters. It is used when we have a fixed set of data and we are looking over the possible values of the parameters (e.g., to fit the model). Once the model is fitted, you now have specific estimates for the parameters so you are generally no longer interested in looking at the likelihood function over the parameters. You can of course substitute your estimated parameters back into the likelihood function; if you estimated using the MLEs then this substitution gives you back the maximised likelihood. Setting that aside, here is how you get the likelihood for a standard linear regression model.

Log-likelihood for a standard linear regression: For a standard linear regression model, with uncorrelated homoscedastic normal error terms, the full model form is:

$$\boldsymbol{Y} = \boldsymbol{x} \boldsymbol{\beta} + \boldsymbol{\varepsilon} \quad \quad \quad \boldsymbol{\varepsilon} \sim \text{N}(\boldsymbol{0}, \sigma^2 \boldsymbol{I}).$$

This model form gives $\boldsymbol{\varepsilon} = \boldsymbol{y} - \boldsymbol{x} \boldsymbol{\beta}$ so that the individual error terms are $\varepsilon_i = y_i - \mathbf{x}_i \cdot \boldsymbol{\beta}$. The error vector represents the deviation of the observed response vector from its expected value under the model. Using this vector gives you the log-likelihood function for the model:

$$\begin{equation} \begin{aligned} \ell_{\boldsymbol{y},\boldsymbol{x}} (\boldsymbol{\beta}, \sigma) &= \sum_{i=1}^n \ln \text{N} ( \varepsilon_i | 0, \sigma^2) \\[6pt] &= \sum_{i=1}^n \ln \text{N} ( y_i - \mathbf{x}_i \cdot \boldsymbol{\beta} | 0, \sigma^2) \\[6pt] &= \frac{n}{2} \ln(2 \pi) + n \ln(\sigma) - \frac{1}{2 \sigma^2} \sum_{i=1}^n (y_i - \mathbf{x}_i \cdot \boldsymbol{\beta})^2 \\[6pt] &= \frac{n}{2} \ln(2 \pi) + n \ln(\sigma) - \frac{1}{2 \sigma^2} || \boldsymbol{y} - \boldsymbol{x} \boldsymbol{\beta} ||^2. \\[6pt] \end{aligned} \end{equation}$$

(The corresponding likelihood function is just the exponential of this.) Maximising this log-likelihood function gives the MLEs $\hat{\boldsymbol{\beta}} = (\boldsymbol{x}^\text{T} \boldsymbol{x})^{-1} (\boldsymbol{x}^\text{T} \boldsymbol{y})$ and $\hat{\sigma}_\text{MLE}^2 = \tfrac{1}{n} || \boldsymbol{y} - \boldsymbol{x} \hat{\boldsymbol{\beta}} ||^2$.$^\dagger$ The coefficient estimator gives the residual vector $\boldsymbol{r} = \boldsymbol{y} - \boldsymbol{x} \hat{\boldsymbol{\beta}} = (\boldsymbol{I} - \boldsymbol{x} (\boldsymbol{x}^\text{T} \boldsymbol{x})^{-1} \boldsymbol{x}^\text{T}) \boldsymbol{y}$. Substitution into the log-likelihood gives the maximised value:

$$\begin{equation} \begin{aligned} \max_{\boldsymbol{\beta}, \sigma} \ell_{\boldsymbol{y},\boldsymbol{x}} (\boldsymbol{\beta}, \sigma) = \ell_{\boldsymbol{y},\boldsymbol{x}} (\hat{\boldsymbol{\beta}}, \hat{\sigma}_\text{MLE}) &= \frac{n}{2} \ln(2 \pi) - n \ln n + n \ln ||\boldsymbol{r}|| - \frac{n}{2}. \\[6pt] &= \frac{n}{2} (\ln(2 \pi)-1) + n \ln \Big( \frac{||\boldsymbol{r}||}{n} \Big). \\[6pt] \end{aligned} \end{equation}$$

You can program this into R or Python if you like, but all of this area of theory has well-known formulae for tests, etc., so you will probably find that you can just look it up in books on regression, rather than using computational optimisers in statistical programming.

---

$^\dagger$ The MLE for the error variance is biased downward. To correct this bias it is standard practice to use the bias-adjusted estimator $\hat{\sigma}^2 = \tfrac{n}{df_{Res}} \cdot \hat{\sigma}^2 = \tfrac{1}{df_{Res}} || \boldsymbol{y} - \boldsymbol{x} \hat{\boldsymbol{\beta}} ||^2$ where $df_{Res}$ is the residual degrees-of-freedom in the regression. This bias-correction is an extension of Bessel's correction for sample variances.