What is the connection between measure theory and statistics?

Measure theory relates to probability theory because "probability" is defined as a normed measure (i.e., a probability measure). In practical statistical work, we don't usually work directly with these measures and the sigma-fields for their domain. Instead we work with density functions that relate to these probability measures. We use the probability measure whenever we go back and compute probabilities for events of interest in the model.

In regression analysis the model form defines a conditional density for the response variable conditional on the explanatory variables, so a typically Gaussian regression model uses the conditional density function:

$$p_{\mathbf{Y}|\mathbf{X}}(\mathbf{y}| \mathbf{x}) = \text{N}(\mathbf{y} | \mathbf{x} \boldsymbol{\beta}, \sigma \mathbf{I}).$$

The resulting probability measure is related to this conditional density by:

$$\mathbb{P}(\mathbf{Y} \in \mathscr{Y} | \mathbf{x}) = \int \limits_\mathscr{Y} p_{\mathbf{Y}|\mathbf{X}}(\mathbf{y}| \mathbf{x}) \ d \mathbf{y} = \int \limits_\mathscr{Y} \text{N}(\mathbf{y} | \mathbf{x} \boldsymbol{\beta}, \sigma \mathbf{I}) \ d \mathbf{y}.$$

Whereas the density function has a single real input argument $\mathbf{y} \in \mathbb{R}^n$ the probability measure has an input argument $\mathscr{Y} \subseteq \mathbb{R}^n$ which is a set of values that is an element of a corresponding sigma-algebra. Typically we are working with real numbers so the sample space is the set $\mathbb{R}^n$ and the sigma-algebra is the corresponding class of Borel sets on $\mathbb{R}^n$. When we are doing statistical analysis we usually don't mention measure theory or the sigma-algebra that the probability measure is defined on, but it is there in the background.

So in a mathematical sense, is the sigma algebra from measure theory related to the variables in $\bf{X}$, or is that just related to the error term $\bf{\epsilon}$?

It may become easier to see if you write the regression as a multivariate normal distributed variable

$$\mathbf{Y} \sim \mathcal{N}(\mathbf{X}\mathbf{B}, \mathbf{\Sigma})$$

which maybe you can see better now, that this is a definition for a probability measure.

A sigma algebra is a set with all intervals/regions $$\prod_{i=1}^n (a_i \leq y_i \leq b_i)$$ and all their their unions. Such intervals describe the events to which we can assign probabilities.

A general, rigorous framework for regression problems can indeed be constructed.

Consider a family of stochastic kernels $\{P_\eta\}_{\eta~\in(a,b)}$ on $(\mathscr Y, \boldsymbol{\mathfrak B}) ,$ that is, for a fixed $B\in\boldsymbol{\mathfrak B}, ~\eta \mapsto P_\eta(B)$ is a measurable map. Let $Y$ be a random variable with values in $\mathscr Y$ and is influenced by regressor $X$ taking values in in $\mathscr X$ with $(\mathscr X, \boldsymbol{\mathfrak A}) .$

Denote by $\mathsf C(\Delta,\mathscr X) $ the class of functions that are continuous in $\theta$ and measurable in $x$ (it is assumed that $\Delta$ is a separable metric space). The regression function $g: \Delta\times \mathscr X\to (a, b), ~g\in \mathsf C(\Delta,\mathscr X)$ determines the distribution of the observation $Y$ given the value of the regressor, that is, if $X=x, $ the distribution of the observation $Y$ is $P_{g_\theta(x) }.$

Regression models involve both cases of regressors: one where $X$ is random and one it being non-random.

When $X$ is random, the conditional distribution of $Y$ given $X=x,$ as discussed above, is given by $\mathsf K_\theta(~\cdot\mid x) =P_{g_\theta(x) }.$ If the marginal distribution of $X$ on $(\mathscr X, \boldsymbol{\mathfrak A}) $ is denoted by $\mu, $ then each independent pair $(Y_i, X_i), ~i\in\overline{1, 2,\ldots, n}, $ has the common joint distribution is $$\mathcal L(Y_i, X_i) =\mathsf K_\theta\otimes \mu.$$

In case of non-random regressors, an experimental design fixes the values $x_i\in\mathscr X;$ the observations $Y_i$s are independent with corresponding distributions being $P_{g_\theta(x_i) }.$

When $\theta, x\in\mathrm{l\! R}^d$ and $g_\theta(x)=\theta^\top x, $ the regression problem obtained is linear. When $g_\theta(x) =\varphi(\theta^\top x), $ where $\varphi$ is continuous, then the regression problem obtained is generalized linear.

--

Reference:

Statistical Decision Theory: Estimation, Testing and Selection, Friedrich Liese, Kalus-J. Miescke, Springer Science$+$Business, $2008, $ sec. $7.5, $ pp. $343-344.$

It's a bit hard to fit in a comment, so I put it as an answer, but we'd be remiss to not mention Gelman's take on this (selected excerpt from his blog):

Stephen Senn quips: “A theoretical statistician knows all about measure theory but has never seen a measurement whereas the actual use of measure theory by the applied statistician is a set of measure zero.”

Which reminds me of Lucien Le Cam’s reply when I asked him once whether he could think of any examples where the distinction between the strong law of large numbers (convergence with probability 1) and the weak law (convergence in probability) made any difference. Le Cam replied, No, he did not know of any examples. Le Cam was the theoretical statistician’s theoretical statistician, so there’s your answer.