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.)
 A: 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)}$$
