Why is the degrees of freedom for a matched pairs $t$-test the number of pairs minus 1?

Cross-posted at Stats SE. It is eligible for bounty in 12 hours, but before I put a bounty on it, I would like to ask it here.

I am used to knowing "degrees of freedom" as $n - r$, where you have the linear model $$\mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\epsilon}$$ with $\mathbf{y} \in \mathbb{R}^n$, $\mathbf{X} \in M_{n \times p}(\mathbb{R})$ the design matrix with rank $r$, $\boldsymbol{\beta} \in \mathbb{R}^p$, $\boldsymbol{\epsilon} \in \mathbb{R}^n$ with $\boldsymbol{\epsilon} \sim \mathcal{N}(\mathbf{0}, \sigma^2 \mathbf{I}_n)$, $\sigma^2 > 0$.

From what I recall of elementary statistics (i.e., pre-linear models with linear algebra), the degrees of freedom for the matched-pairs $t$-test is the number of differences minus $1$. So this would entail $\mathbf{X}$ having rank 1, perhaps. Is this correct? If so, what would the $\mathbf{y}=\mathbf{X}\boldsymbol{\beta}+\boldsymbol{\epsilon}$ equation be? If not, why is $n-1$ the degrees of freedom for the matched-pairs $t$-test?

To understand the context, suppose I have a mixed-effects model $$y_{ijk} = \mu_i + \text{ some random effects} + e_{ijk}$$ where $i = 1, 2$, $j = 1, \dots, 8$, and $k = 1, 2$. There is nothing special about $\mu_i$ other than that it's a fixed effect, and $e_{ijk} \overset{iid}{\sim}\mathcal{N}(0, \sigma^2_e)$. I'm assuming that the random effects are irrelevant to this problem, since we only care about the fixed effects in this case. (For those of you not familiar with what random effects are, they are just a bunch of independent normally-distributed random variables with mean 0 and their own respective variance. Fixed effects are just constants.)

I would like to provide a confidence interval for $\mu_1 - \mu_2$.

I have already shown that $\bar{d}_\cdot = \dfrac{1}{8}\sum d_j$ is an unbiased estimator of $\mu_1 - \mu_2$, where $d_j = \bar{y}_{1j\cdot} - \bar{y}_{2j\cdot}$, $\bar{y}_{1j\cdot} = \dfrac{1}{2}\sum_{k}y_{1jk}$, and $\bar{y}_{21\cdot}$ is defined similarly. The point estimate $\bar{d}_{\cdot}$ has been computed.

I have already shown that $$s^2_d = \dfrac{\sum_{j}(d_j - \bar{d}_{\cdot})^2}{8-1}$$ is an unbiased estimator of the variance of $d_j$, and thus, $$\sqrt{\dfrac{s^2_d}{8}}$$ is the standard error of $\bar{d}_{\cdot}$. This has been computed.

Now the last part is figuring out the degrees of freedom. For this step, I usually try to find the design matrix - which obviously has rank 2 - but I have the solution to this problem, and it says that the degrees of freedom is $8-1$.

In the context of finding the rank of a design matrix, why are the degrees of freedom $8-1$?

Edited to add: Perhaps helpful in this discussion is how the test statistic is defined. Suppose I have a parameter vector $\boldsymbol{\beta}$. In this case, $$\boldsymbol{\beta} = \begin{bmatrix} \mu_1 \\ \mu_2 \end{bmatrix}$$ (unless I'm missing something entirely). We are essentially performing the hypothesis test $$\mathbf{c}^{\prime}\boldsymbol{\beta} = 0$$ where $\mathbf{c}^{\prime} = \begin{bmatrix} 1 & -1 \end{bmatrix}$. Then, the test statistic is given by $$t = \dfrac{c^{\prime}\hat{\boldsymbol{\beta}}}{\sqrt{\hat{\sigma}^2c^{\prime}(\mathbf{X}^{\prime}\mathbf{X})^{-1}\mathbf{c}}}$$ which would be tested against a central $t$-distribution with $n - r$ degrees of freedom, where $\mathbf{X}$ is the design matrix as above, and $$\hat{\sigma}^2 = \dfrac{\mathbf{y}^{\prime}(\mathbf{I}-\mathbf{P}_{\mathbf{X}})\mathbf{y}}{n-r}$$ where $\mathbf{P}_{\mathbf{X}} = \mathbf{X}(\mathbf{X}^{\prime}\mathbf{X})^{-1}\mathbf{X}^{\prime}$. This is just how I know it to be defined from my coursework; I may have this wrong. However, in general, $\boldsymbol{\beta}$ is the vector from the definition of the linear model equation above.

The idea is this: let $$\mathbf{y} = \begin{bmatrix} d_1 \\ \vdots \\ d_n \end{bmatrix}$$ and $\boldsymbol{\beta} = [\mu_1 - \mu_2]$. Then our linear model is then $$\mathbf{y} = \mathbf{1}_{n \times 1}\boldsymbol{\beta} + \boldsymbol{\epsilon}$$ where $\mathbf{1}_{n \times 1}$ is the $n$-vector of all ones, and $$\boldsymbol{\epsilon} = \begin{bmatrix} \epsilon_1 \\ \vdots \\ \epsilon_n \end{bmatrix} \sim \mathcal{N}(\mathbf{0}, \sigma^2\mathbf{I}_n)\text{.}$$ Obviously $\mathbf{X} = \mathbf{1}_{n \times 1}$ has rank $1$, so then we have $n-1$ degrees of freedom.
How do we know to set $\boldsymbol{\beta}$ equal to $[\mu_1 - \mu_2]$? Recall that $$\mathbb{E}[\mathbf{y}] = \mathbf{X}\boldsymbol{\beta}$$ and as it can be easily seen, $\mathbb{E}[d_j] = \mu_1 - \mu_2$ for all $j$. Given our $\mathbf{X}$, it is obvious what $\boldsymbol{\beta}$ should be.
Set $\mathbf{c}^{\prime} = [1]$. Then our hypothesis test is $$H_0: \mathbf{c}^{\prime}\boldsymbol{\beta} = 0\text{.}$$ Our test statistic is thus $$\dfrac{\mathbf{c}^{\prime}\hat{\boldsymbol{\beta}}}{\sqrt{\hat{\sigma}^2\mathbf{c}^{\prime}\left(\mathbf{X}^{\prime}\mathbf{X}\right)^{-1}\mathbf{c}}}\text{.}$$ We have $$\hat{\sigma}^2 = \dfrac{\mathbf{y}^{\prime}(\mathbf{I}-\mathbf{P}_{\mathbf{X}})\mathbf{y}}{n-r(\mathbf{X})}\text{.}$$ After some work, it can be shown that $$\mathbf{P}_\mathbf{X} = \mathbf{P}_{\mathbf{1}_{n \times 1}} = \mathbf{1}_{n \times 1}\left(\dfrac{1}{n}\right)\mathbf{1}^{\prime}\text{.}$$ It can also be shown that $\mathbf{I}-\mathbf{P}_{\mathbf{X}}$ is symmetric and idempotent. So, $$\hat{\sigma}^2 = \dfrac{\left\|\left[\mathbf{I}-\mathbf{1}_{n \times 1}\left(\dfrac{1}{n}\right)\mathbf{1}^{\prime}\right]\mathbf{y} \right\|^2}{n-1} = \dfrac{\left\|\begin{bmatrix} d_1 \\ \vdots \\ d_n \end{bmatrix} - \begin{bmatrix} \bar{d}_\cdot \\ \vdots \\ \bar{d}_\cdot \end{bmatrix} \right\|^2}{n-1} = \dfrac{\sum_{i=1}^{n}(d_i-\bar{d}_{\cdot})^2}{n-1} = s^2_d$$ and $$\mathbf{X}^{\prime}\mathbf{X} = \mathbf{1}_{n \times 1}^{\prime}\mathbf{1}_{n \times 1} = n$$ which obviously has inverse $1/n$, thus giving a test statistic $$\dfrac{\hat\mu_1-\hat\mu_2}{\sqrt{s^2_d/n}}$$ which would be tested on a central $t$-distribution with $n - 1$ degrees of freedom as desired.