The uncertainity about the weight matrix in linear model

Question

Textbook literature often denotes the estimated weight matrix of a linear regression model $y = Wx + \epsilon, \epsilon \sim \mathcal{N}(0,\sigma^2)$ by $\hat W$ due to the inherent variability in the sample size. Limited access to all possible instances hinders our ability to find the true matrix $W$. Additionally, the relationship between the inputs ($x$) and the output ($y$) might not be linear, resulting in uncertainty represented by the random variable $\epsilon$. Incorporating both sources of randomness into the estimate of $W$, we can see that $\epsilon$ is random $\rightarrow y$ is random $\rightarrow \hat W$ is random. Therefore, $\hat W | data \sim \mathcal{N}(W,\sigma^2(X^TX)^{-1})$, where $X \in \mathbb{R}^{n \times p}$, and $p$ is the number of features while $n$ is the sample size.

1. Do you agree with me that in this model, we have two sources of randomness/uncertainty: the sample size variation and the uncertainty induced by the input-output relationship?
2. As $n \rightarrow \infty$, $(X^TX)^{-1}$ approaches the true population inverse covariance matrix of the feature vector, if the features are naturally random. However, if the features are deterministic (hence the randomness of the input is due to the sample size, not to inherent randomness in the input), then I guess each element of the covariance matrix of the estimate of $W$ should go to zero, indicating that the uncertainty about the model weight matrix becomes negligible. However, I think it will converge to constant values, but what do those values mean? In this case, I think the model projects spurious randomness, and this should be appropriately treated.
3. If the distribution of $\hat W$ accounts for all sources of randomness, Why do we need a Bayesian approach to estimate the weight then? There is already a way to estimate the uncertainty when predicting a new instance: $p(y_{test} | Data) = \int p(y_{test} | \hat W,x_{test})p(\hat W|data) d \hat W$, where $p(\hat W|data)\sim \mathcal{N}(W,\sigma^2(X^TX)^{-1}$) (assuming that I do have the model $p(y_{test} | \hat W,x_{test})$. I understand that the original assumption in this modelling is that we have true value of $W$ which is the true mean of the distribution of $\hat W$ and based on this assumption, we should only use that value. However, I am a bit confused of why we can't model the uncertainty this way.

Answer 1

1. Notice that the sample covariance is $S_n = \frac{1}{n}\sum_i x_i x_i^T=\frac{1}{n}X^TX$. If the features are drawn from a distribution with covariance matrix $\Sigma$, then $S_n$ converges in probability to $\Sigma$ when $n$ approaches infinity. In other words, $(X^TX)^{-1} \to n^{-1} \Sigma^{-1} \to 0$ when $n \to \infty$.

2. In linear regression the features are usually considered as fixed. In particular, the distribution of the estimator $\hat W \sim \mathcal N(W,\sigma^2(X^TX)^{-1})$ is due to the randomness of $y$ only. (Intuitively, it is the distribution you will get by sampling $y$ over and over with the same features). If you would like to consider how $\hat W$ varies due to randomness of the features, you will have to additionally model the distribution of $x$, which will result in a much more complicated distribution of $\hat W$ (and probably impossible to find in closed form).

3. The estimator $\hat W$ is fixed by the data, so the expression $p(\hat W|data)$ doesn't make a lot of sense (it is just equal to 1). To make sense of it you would have to treat $W$ itself as the random variable, which is exactly the Bayesian approach.