Where does the t-distribution come from to calculate the p-value in linear regression? I apologize in advance if I'm fundamentally misunderstanding something here.  But I'm trying to figure out where the data is produced that forms the t-distribution on which one calculates the p-value for a coefficient in linear regression.  Is the regression run multiple times, maybe on bootstrapped samples, in order to create this distribution?
 A: Under the assumptions of the linear regression model, i.e., that
$$Y_i = \beta_0 + \beta_1X_{i1} + \cdots + \beta_pX_{ip} + \epsilon_i,$$
where $\epsilon_i \sim \text{N}(0, \sigma^2),$ the estimator of the regression coefficients $\boldsymbol{\hat\beta}$ is by construction a linear combination of normally distributed random variables, and is therefore also normally distributed. To test the assumption that a given regression coefficient $\beta_j$ is zero, we use the test statistic
$$T = \frac{\hat \beta_j}{\hat\sigma_{\hat\beta_j}}.$$
If we pretend that $\hat \sigma_{\hat\beta_j}$ is the true standard deviation of $\hat \beta_j$, we get an approximate test where we compare the value of $T$ to the quantiles of the standard normal distribution. The problem is that $\epsilon_1, \cdots, \epsilon_p$ factors into the estimate of $\sigma_{\hat\beta_j}.$ It turns out that the distribution of the estimator $\hat\sigma_{\hat\beta_j}$ has a
$\chi^2$- distribution. The ratio of a normally distributed random variable and a $\chi^2$-distributed random variable follows a $t$-distribution. The regression is not performed multiple times, the $p$-values come from the quantiles of the $t$-distribution because, under the model, $T$ follows a $t-$distribution. There is no additional data, $T$ is a function of the observed data.
A: Suppose $Y_i = \alpha + \beta x_i +\varepsilon_i$ for $i=1,\ldots,n$ and $\alpha,\beta, x_i$ are not random and $\varepsilon_1,\ldots,\varepsilon_n \sim \text{i.i.d.} \operatorname N(0,\sigma^2).$ The least-squares estimators $\widehat\alpha, \widehat\beta$ of $\alpha,\beta$ are given by
\begin{align}
\widehat\beta & = \frac{\sum_{i=1}^n (x_i-\overline x)(Y_i-\overline Y)}{\sum_{i=1}^n (x_i-\overline x)^2} \\[10pt]
\widehat\alpha & = \overline Y - \widehat\beta\,\overline x \\[10pt]
\text{where } \overline x & = \frac{\sum_{i=1}^n x_i} n \text{ and } \overline Y = \frac{\sum_{i=1}^n Y_i} n.
\end{align}
Now observe that $\widehat\beta$ is a linear combination of $Y_1,\ldots,Y_n$ and the coefficients in the linear combination are not random. Thus $\widehat\beta$ has a normal distribution. By using those coefficients and doing some algebra, you see that
$$
\widehat\beta \sim \operatorname N\left( \beta, \frac{\sigma^2}{\sum_{i=1}^n (x_i-\overline x)^2} \right).
$$
The residuals are
$$
\widehat\varepsilon_i = Y_i - \widehat Y_i = Y_i - (\widehat\alpha + \widehat\beta x_i).
$$
This is also a linear combination of $Y_i,\,\, i=1,\ldots,n,$ so that this is normally distributed, and it should be clear that $\operatorname E\widehat\varepsilon_i=0.$ But $\operatorname{var}(\widehat\varepsilon_i)$ and $\operatorname{cov}(\varepsilon_j,\varepsilon_k)$ depend on $x_i\,\, i=1,\ldots, n$ in a seemingly complicated way. But some routine if laborious algebra will show you that $\operatorname{cov}\left(\,\widehat\varepsilon_i, \widehat\beta\,\right) = 0.$ Given the initial assumptions about distributions, this implies that the vector $\left(\widehat\varepsilon_i : i=1,\ldots,n\right)$ and the vector $\left(\widehat\alpha,\widehat\beta\right)$ are independent of each other.
This next step I don't know off-hand how to do without talking about geometry. Observe that


*

*the vector $\left( \widehat\varepsilon_i: i=1,\ldots,n\right)$ is the orthogonal projection of the vector $\left( Y_i: i=1,\ldots,n\right)$ onto the orthogonal complement of the space spanned by the two vectors $\big(\, \underbrace{\,1,\ldots,1\,}_{n \text{ components}}\,\big),$ $\big(x_1,\ldots,x_n\big).$ And

*This orthogonal projection maps $\operatorname E(Y_i : i=1,\ldots, n)$ to the zero vector.
A consequence is that $\left( \widehat\varepsilon_i: i=1,\ldots,n\right)$ has a spherically symmetric normal distribution with expected value $0$ in an $(n-2)$-dimensional space.
Therefore
$$
\frac 1 {\sigma^2} \sum_{i=1}^n \widehat\varepsilon_i^2 \sim \chi^2_{n-2}.
$$
So we have
$$
\frac{(\widehat\beta-\beta)\left/\left(\sigma\left/\sqrt{\sum_{i=1}^n (x_i-\overline x)^2} \right. \right)\right.}{(1/\sigma)\sqrt{\sum_{i=1}^n \widehat\varepsilon_i^2/(n-2)}}
= \frac{(\widehat\beta-\beta)\sqrt{\sum_{i=1}^n (x_i-\overline x)^2}}{\sqrt{\sum_{i=1}^n \widehat\varepsilon_i^2/(n-2)}} = \frac{Z}{\sqrt{\chi^2_{n-2}/(n-2)}}
$$
and the numerator and denominator are independent. Thus this has a t-distribution, and can accordingly be used to find a confidence interval for $\beta.$
