# How to calculate p-value for multivariate linear regression

Software packages that calculate regressions sometimes also return p-values. I want to understand how to calculate this p-value by hand.

Here's what I think I understand:

I want to calculate the probability that ${\bf y} = {\bf X\beta}$ where $\bf y$ is a column vector of $m$ dependent-variable values, $\bf X$ is an $m\times n$ matrix where each row is a set of observations of $n$ variables with normally distributed errors, and $\bf\beta$ is a column vector of length $m$.

In particular to calculate the p-value I assume that the regression coefficients $\bf\beta$ are zero. So I want to calculate the probability that ${\bf y} = \bf \bf\epsilon$, where the $\bf\epsilon$ is a vector of normally distributed errors with means = 0? Is this correct?

If so, I am not sure what the next step is for calculating the p-value. (Among other things, where do I get the variance parameters for the normal distributions?)

(Feel free to point me to an answer to an existing question, obviously. I haven't yet found an existing question with an answer to mine, but I won't be surprised if one exists.)

• You can have a p-value for an F-test, comparing the differences in residuals with and without the variable (anova), or for a t-test, calculating/estimating the variance for the parameters and use this to calculate t-scores. – Martijn Weterings Jun 20 '18 at 18:24
• Roughly, yes. It will be a lot easier, however, if you think it through in the case where there is only one $\beta$ (and no intercept). Suppose you are trying to predict weight from height, and we know zero height is zero weight. Intuitively, once you do the regression, you can 'predict' what the weight 'should' be for each height in your data. This will be different, in most cases, from the actual data for weight you started with. So, you can estimate the variance of the errors (which, in theory, is what makes the actual and predicted be different). – eSurfsnake Jun 20 '18 at 19:41
• Are you looking for one p-value for the model as a whole? – Sal Mangiafico Jun 20 '18 at 23:36
• @SalMangiafico I was thinking of the overall p-value, but on reflection I realize that I actually want to understand how to calculate the per-variable p-values. – Mars Jun 21 '18 at 2:49
• @MartijnWeterings I was unaware of the usual role of F-tests and t-tests in significance testing. I think I need to read up on them. As you can see, I am quite ignorant about statistics. (I'm comfortable with probability, though.) – Mars Jun 21 '18 at 2:51

## t-test

With a t-test you standardize the measured parameters by dividing by them by the variance. If the variance is an estimate then this standardized value will be distributed according to the t-distribution (otherwise, if the variance of the distribution of the errors is known, then you have a z-distribution)

Say your measurement is:

$$y_{obs} = X\beta + \epsilon \quad \text{with} \quad \epsilon \sim N(0,\sigma^2*I)$$

Then your estimate $\hat\beta$ is:

$$\begin{array}\\ \hat\beta & = & (X^TX)^{-1}X^T y_{obs} \\ & = &(X^TX)^{-1}X^T (X\beta + \epsilon) \\ & = & \beta + (X^TX)^{-1}X^T \epsilon \end{array}$$

So your estimate $\hat\beta$ will be the true vector $\beta$ plus a term based on the error $\epsilon$. If $\epsilon \sim N(0,\sigma^2I)$ then

$$\hat\beta \sim N(\beta,(X^tX)^{-1}\sigma^2)$$

Note: I can not make the change of the $(X^TX)^{-1}X$ term into $(X^TX)^{-1}$ intuitive, but to derive this you would express $\text{Var}(\hat\beta) = \text{Var}((X^TX)^{-1}X^T\epsilon) = (X^TX)^{-1}X^T \, \sigma^2I \, ((X^TX)^{-1}X^T)^T$ and eliminate some of those terms

The unknown $\sigma$ will be estimated based on the sum of squares of the residuals multiplied by the ratio of the degrees of freedom in the residual terms and the total number of measurements/error-terms (in a similar fashion as Bessel's correction in the corrected sample variance) .

Then from this point you can pick up the expression of p-values for single $H_j: \beta_j = 0$ as standard t-tests (although due to the possible correlation in the distribution of different $\beta_j$, different more powerfull tests than individual t-tests could be done).

## F-test

With the F-test you use the F-distribution which describes the ratio of two chi-squared distributed variables. This works as a hypothesis test when we compare the variance of a model and residuals (both are chi-square distributed when we assume that a certain model parameter $\beta_j$ has no effect)

The residual term of a model has $n-p$ degrees of freedom, with $n$ the number of observations/errors and $p$ the number of parameters that are used to fit the model. You could see this intuitively as the residuals being obtained from the errors by projecting the errors onto the space from the columns perpendicular to the model columns $X$ (this space has dimension $n-p$). A projection of a multivariate normal distributed variable is itself a multivariate normal distributed variable, but with a lower dimension. So while you may have $n$ residuals. They are actually $n-p$ residuals embedded in an $n$ dimensional space.

Now when you consider adding an extra variable to model 1, to obtain model 2, then you could analyze this by considering the errors being projected onto a smaller space. If the model 2 has no effect (ie. the added columns to make model $X_2$ from $X_1$ are just random) then one could state a null hypothesis that the reduced sum of squared residuals for model 1 and model 2 are equal. This is what is tested in an F-test (using the ratio of those reduced residuals) to obtain a p-value for the effect of changing model 1 into 2 (and you could do this for every variable $\beta_i$ where the way that you do this changes a bit see for instance How to interpret type I, type II, and type III ANOVA and MANOVA?).

So you split the sum of squared residuals of a simple model $RSS_{simple}$ into two projections (representing independent variables if the null hypothesis is true). One part is a projection onto the (smaller) space of a full model $RSS_{full}$ and the other part is the projection onto the space spanned by the model (which can be expressed by the difference) $RSS_{simple}-RSS_{full}$. And the ratio used in the F-test is

$$F = \frac{\left(\frac{RSS_{simple}-RSS_{full}}{p_{full}-p_{simple}}\right)}{\left(\frac{RSS_{full}}{n-p_{full}}\right)}$$

• Thanks Martijn! The t-test explanation was very clear and thorough. There's a lot I don't understand in the F-test section, but I don't blame your explanation for that. I need to do some background reading, and I don't think an answer here has to fill in all gaps. – Mars Jun 21 '18 at 16:12
• A few questions about the F-test section: Model 1 is the one in which we assume $\beta_j$ has no effect, and model 2 is the model in which we assume it does, right? The projection is onto a smaller space in this sense (?): The model 1 space has $k$ dimensions, and by assuming a specific value for $\beta_j$ we constrain values, and project them onto, a subspace defined by the value of $\beta_j$. (I was confused at first because I thought that adding a parameter added a dimension.) – Mars Jun 21 '18 at 16:17
• The projection is onto a smaller space is a bit ambiguous. You can consider the space spanned by the column vectors of $X$. In this space you have for instance the vector $\hat y$. And the space that is composed of the vectors perpendicular to this (the orthogonal complement). In this space you have for instance the residual vector $y- \hat y$. For the smaller model the 'model space' is smaller but the 'residual space' is larger. – Martijn Weterings Jun 21 '18 at 17:32